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EXECUTIVE  SUMMARY 


The  development  of  a  numerical  model  of  steady-state  plow/soil  interaction  in  support  of 
cable  burial  operations  on  the  seafloor  was  addressed.  It  was  guided  by  previous  laboratory 
experiments  showing  a  linear  relationship  between  plow  resistance  and  steady-state  velocity  up 
to  a  critical  velocity  in  the  range  of  10  fps,  beyond  which  resistance  tended  to  remain  constant. 
Further,  the  data  showed  that  soil  stress  and  pore  water  flow  near  the  plow  were  three- 
dimensional  quantities.  Nonetheless,  because  of  the  unaffordable  expense  of  full  three- 
dimensional  simulations,  two-dimensional  simulations  were  used  to  demonstrate  behavior  in  a 
representative  horizontal  plane.  Coupled  field  equations  for  equilibrium  of  an  elasto-plastic 
medium  and  pore  water  diffusion  were  solved  in  the  presence  of  rigid  surfaces  contacting  and 
sliding  through  fully  saturated  soil.  Commercial  finite  element  technology  imposed  a  multistep 
transient  solution  approach,  and  the  representative  software  employed  in  the  solution  further 
imposed  the  usage  of  high  order  finite  element  formulations  for  soil.  These  restrictions 
contributed  to  a  computationally  intensive  model,  prohibiting  a  three-dimensional  solution  using 
available  engineering  workstation  computers. 

A  summary  was  also  provided  of  two  alternative  approaches:  (1)  a  simplified  semi- 
empirical  model,  and  (2)  a  new  single-step  finite  element  demonstration  model  for  solving  the 
steady  state  problem  directly.  The  former  model  included  three-dimensional  effects.  The  latter 
model  exploited  equivalency  of  time  and  spatial  derivatives  in  the  formulation  of  the  problem. 
Both  models  were  successful,  but  the  latter  demonstrated  potential  for  ultimately  achieving  a 
rational  numerical  model  in  three  dimensions. 

The  multistep  transient  solution  considered  a  simplified  two-dimensional  plane  strain 
formulation  of  a  rigid  thin-blade  plow  interacting  with  soil.  Though  inherently  incapable  of 
considering  total  plow  resistance  force,  the  model  featured  finite  strains,  a  modified  Drucker- 
Prager  plasticity  constitutive  model,  and  large  sliding  along  a  frictional  plow/soil  interface. 
Special  handling  techniques  enabled  reasonable  responses  for  variables  and  workable  compute 
times. 

Time  histories  of  unit  plow  forces,  in  lieu  of  total  plow  forces,  displayed  distinct  transient 
and  steady-state  phases  reminiscent  of  total  force  histories  recorded  in  the  laboratory,  and  the 
computed  relationship  between  unit  plow  force  and  velocity  remained  linear  from  zero  to  13  fps, 
although  no  critical  velocity  was  evident  in  the  results.  Permanent  soil  deformation  in  the  wake 
of  the  plow  replicated  laboratory  observation.  Concentrations  of  effective  stress  and  contact 
forces  distributed  over  the  plow  blade  were  consistent  with  patterns  of  severe  wear  on  plow 
blades  observed  in  the  field.  Computed  soil  dilatation  response  was  reasonable,  and  negative 
pore  pressures  were  representative  of  measured  response,  although  a  provision  for  modeling 
cavitation  in  the  pore  fluid  did  not  present  itself.  The  model  verified  that  only  the  tip  of  the  thin- 
blade  plow  contacted  the  soil,  it  complemented  available  experimental  data  regarding  lateral 
plow  force,  and  it  also  evidenced  expansive  soil  behavior  as  should  be  anticipated  for  dilatant 
seafloor  soil .  Moreover,  the  model  explained  how  the  mean  pressure  stress  builds  forward  of  the 
plow  tip  strengthening  the  pressure-sensitive  soil  through  which  it  penetrates,  and  suggests  why 
past  field  measurements  of  resistance  exceeded  predictions  from  conventional  analysis 
procedures. 
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A  parameter  study  further  demonstrated  the  engineering  utility  of  rational  numerical 
plow/soil  interaction  models.  The  consequences  of  coulomb  friction  at  the  plow/soil  interface 
were  potentially  significant  for  tow  forces  but  not  for  lateral  forces,  and  while  the  potential  for 
blade  sharpness  to  reduce  tow  force  was  also  significant,  sharper  blades  also  reduce  permanent 
soil  displacement  and  the  potential  effectiveness  of  trench  formation. 

Commercial  nonlinear  finite  element  technology  demonstrated  substantial  capability  and 
surprising  resilience  in  a  highly  specialized  application  to  plow/soil  interaction  modeling.  Future 
development  of  affordable  three-dimensional  numerical  models  should  also  consider  special- 
purpose  finite  element  software  solutions  based  on  a  steady-state,  single-step  finite  element 
method  so  that  solutions  can  be  realized  for  engineering  design  applications. 
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1.0  INTRODUCTION 


The  problem  addressed  in  this  report  is  the  interaction  of  a  rigid  thin-blade  plow  with 
dense  100  percent  saturated  sand  as  the  plow  is  towed  at  a  constant  rate  along  the  surface  of  the 
seafloor  for  the  purpose  of  creating  a  trench.  The  important  engineering  information  sought  is 
the  resisting  force  acting  on  the  blade  as  it  cuts  through  the  soil.  Plow  resistance  force  or  tow 
force  and  energy  consumption  measured  in  the  field  have  reportedly  far  exceeded  values 
predicted  by  conventional  theories  and  by  previous  laboratory  model  tests.  Reasons  for  the 
discrepancies  include  increase  in  soil  strength  related  to  plow  velocity,  and  in  laboratory 
investigations,  inadequate  preparation  of  soil  test  conditions  which  were  not  representative  of 
fully  saturated  seafloor  soils. 

A  substantial  laboratory  experimental  program  was  developed  at  the  Naval  Facilities 
Engineering  Service  Center  (NFESC)  to  resolve  the  laboratory  issues.  The  laboratory  apparatus 
for  measuring  resistance  forces  on  scale  models  of  thin  blade  plows  is  depicted  in  Figure  1-1. 
This  sketch  provides  an  overview  of  the  physical  problem  with  which  the  engineering 
investigation  into  the  mechanics  of  plow/soil  interaction  is  concerned.  The  portable  soil  test 
fixture  shown  in  Figure  1 -la  has  a  controlled  hydraulic  ram-cable  linkage  system  which  pulls 
with  a  force  sufficient  to  maintain  a  preset  constant  velocity  for  the  trolley  and  test  article  or 
thin  blade  plow.  The  instrument  package  measures  the  resistance  force,  Fr,  as  well  as  the 
moment  ofFR.  In  addition  to  providing  data  at  one  atmosphere,  the  apparatus  is  also 
transportable  for  testing  at  various  simulated  ocean  depths  at  distant  hyperbaric  chamber  test 
sites. 

The  geometry  of  a  typical  thin  blade  plow  is  shown  in  Figure  1-lb.  The  main 
characteristics  are  the  blade’s  thickness,  length,  depth,  sweep  angle,  and  tip  included  angle.  For 
the  plow  shown  these  values  are,  respectively,  1/2  inch,  6.5  inches,  9-inch  nominal,  90  and  60 
degrees. 

The  experimental  program  has  produced  two  main  technical  reports  to  date  on  test 
results  with  two  types  of  fully  saturated  sand  materials  —  Cable  et  al.  (1993)  and  Girard  and 
Taylor  (1994).  Typical  primary  data  obtained  from  these  tests  are  shown  in  Figure  l-2a,  for 
6-inch-deep  thin  blade  plow  test  articles  having  various  thicknesses,  and  under  various 
conditions  of  soil  saturation  (the  full  100  percent  saturation  condition  is  of  primary  interest).  In 
the  range  of  test  velocities  up  to  about  8  fps,  it  is  seen  that  the  relationship  between  plow 
resistance  and  plow  velocity  is  linear.  The  soil  type  was  a  medium  dense  grain-size  sand  in 
these  data,  i.e.,  a  dense  Yorba  River  sand.  Similar  results,  in  Figure  l-2b,  for  the  9-inch-deep 
thin  blade  plow  shows  the  relationship  to  be  linear  up  to  8  fps  with  a  steeper  slope.  Beyond 
these  so-called  critical  velocities,  the  behavior  tends  toward  a  constant  resistance. 

To  help  understand  the  controlling  physics  of  this  problem,  the  development  of  a 
numerical  model  of  plow/soil  interaction  was  authorized.  How  well  commercial,  general 
purpose  finite  element  technology  would  perform  in  providing  this  capability  was  an  immediate 
goal  of  the  study.  The  aforementioned  experimental  program  was  instrumental  in  guiding  the 
development. 

Numerical  modeling  of  the  coupled  equilibrium  and  porous  media  flow  fields  in  the 
presence  of  rigid  surfaces  and  sliding  contact  is  a  complex  simulation  problem.  However, 
solutions  are  feasible  although  they  generally  cannot  be  rushed.  It  is  mandatory  to  maintain 
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control  of  the  development  each  step  of  the  way  by  knowing  the  theoretical  and  numerical 
meaning  of  model  data.  In  this  way,  a  reliable  and  rational  numerical  model  can  provide  new 
data  and  detailed  insight  describing  the  behavior  of  plow/soil  interaction. 

Development  of  a  numerical  model  of  plow/soil  interaction  using  available  commercial 
finite  element  technology  was  pursued.  The  approach  taken  includes  modeling  both  nonlinear 
saturated  soil  material  behavior  and  the  mechanics  of  plow/soil  interaction. 


1.1  Objectives 

The  overall  objective  is  to  provide  a  numerical  model  for  plow/soil  interaction,  and  to: 

1 .  Predict  the  resistance  of  a  thin  blade  plow  to  steady-state  motion  through  seafloor 
soil. 

2.  Provide  an  analysis  tool  for  evaluating  system  parameters  on  a  rational  basis. 

3.  Assess  the  capability  of  commercial,  general-purpose  finite  element  computer 
programs  to  provide  these  capabilities. 


1.2  Scope 

ABAQUS/Standard,  Version  5.3  (Hibbitt,  Karlsson,  and  Sorensen,  Inc.,  (HKS),  1993), 
was  the  computer  program  employed  while  conducting  the  study.  This  program  is  well-regarded 
for  solving  nonlinear  problems  in  structural  and  solid  mechanics  and  is  very  representative  of 
commercial  technology.  With  regard  to  geotechnical  problems,  it  features  special  capabilities  for 
effective  stress  analysis  and  three  alternative  plasticity-based  soil  models: 

1 .  Modified  Drucker-Prager  Model 

2.  Modified  Cap  Model 

3 .  Modified  Cam-Clay  Model 

This  program  also  has  substantial  interface  modeling  features  including  sliding  and  friction 

which  are  also  applicable  to  geotechnical  problems. 

ABAQUS/Standard  is  an  implicit  finite  element  computer  program  suitable  for  quasi¬ 
static  or  dynamic  problems.  Like  other  commercial  finite  element  programs,  it  has  no  capability 
for  solving  steady-state  problems  directly.  That  is,  the  user  must  simulate  a  steady-state  problem 
using  a  multistep  transient  solution  approach.  The  steady-state  problem  is  viewed  as  a  special 
case  of  the  more  general  problem,  requiring  the  solution  of  an  initial  transient  phase  preceding 
the  solution  of  the  steady-state  phase. 

Appendix  A  contains  a  discussion  of  some  relevant  technical  issues  and  the  consequences 
of  project  constraints  which  were  prepared  early  on  in  this  study. 
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Figure  1-1.  Laboratory  testing  apparatus. 


Plow  Resistance  vs.  Plow  Velocity 
Comparison  of  6"  Blades 
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Figure  1-2.  Plow  resistance  vs.  plow  velocity. 


Plow  Resistance  vs.  Plow  Velocity 
9"  Blade  Comparison 
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(b)  Comparison  of  9-inch  blades. 


1.3  Apparent  Technical  Risk 

No  general  purpose  nonlinear  finite  element  computer  program  is  likely  to  have  all  the 
technical  capabilities  desired  to  solve  a  particularly  complex  problem.  In  most  cases,  solutions 
are  coerced,  at  least  to  some  degree,  from  the  program  selected.  Often  technical  difficulty  simply 
cannot  be  avoided  because  solution  requirements  simply  exceed  the  resources  of  commercial 
state-of-the-art  computational  mechanics.  Four  areas  of  numerical  modeling  complexity 
encountered  in  the  present  treatment  of  coupled  equilibrium/porous  media  flow  problems  are 
discussed  in  1.3.1,  1.3.2,  1.3.3,  and  1.3.4.  They  are  generally  related  to  the  use  of  high  order 
elements  and  constitutive  models. 

1.3.1  High  Order  Element  Burden.  Effective  stress  soil  models  are  required  so  that 
pore  water  pressure  effects  may  be  included  in  the  analysis.  The  effective  stress  principle  in  soil 
mechanics  provides  for  this  by  conceptualizing  the  mechanical  behavior  of  the  soil  skeleton  in 
terms  of  an  effective  stress,  while  the  overall  equilibrium  equations  for  the  soil  are  formulated  in 
terms  of  the  total  stress  (algebraic  sum  of  effective  stress  and  pore  pressure).  In  this  way,  pore 
pressure  may  be  interpolated  independently  from  displacements  in  a  corresponding  finite  element 
spatial  discretization  of  the  problem.  Most  effective  stress  finite  element  formulations  used 
throughout  the  industry  work  in  this  way  (Zienkiewicz  and  Taylor,  1989)  i.e.,  interpolation 
functions  for  the  pore  pressure  field  and  the  soil  skeleton  displacement  field  are  different, 
whether  the  soil  element  is  regarded  as  a  high  order  or  low  order  formulation. 

Usage  of  a  high  order  instead  of  a  low  order  finite  element  formulation  for  analysis  of 
coupled  equilibrium/porous  media  flow  was  imposed  by  the  software  employed  in  this  study. 
This  lead  to  lengthy  compute  times  and  large  input  and  output  data  files.  Bearing  in  mind  that 
both  displacements  and  pore  pressures  are  evaluated  at  comer  nodes  in  this  formulation,  the  two- 
dimensional  version  has  four  midside  and  four  comer  nodes  and  the  three-dimensional  version 
has  twelve  midside  and  eight  comer  nodes.  Only  displacements  are  evaluated  at  midside  nodes. 
This  is  in  comparison  with  four  and  eight  comer  nodes,  respectively,  for  two-  and  three- 
dimensional  elements  of  low  order  finite  element  formulations.  Low  order  elements  have  no 
midside  nodes,  though  there  is  one  center  node  for  evaluating  pore  pressure.  Further,  the  high 
order  two-dimensional  element  has  9  internal  points  (Gauss  points)  where  stress  and  strain  are 
evaluated,  and  also  4  points  where  velocity  gradient  is  evaluated.  The  corresponding  low  order 
element  has  only  4  internal  points  for  stress/strain  evaluation  and  only  a  single  point  for 
evaluating  velocity  gradient.  A  low  order  soil  element  formulation  would  have  been  a  welcome 
alternative  in  this  study. 

1.3.2  Clarity  of  Boundary  Condition  Specification.  The  element  formulation  for 
coupled  equilibrium/porous  media  flow  is  necessarily  a  mixed  finite  element  formulation;  as 
indicated  above,  displacements  and  pore  pressures  are  the  primary  variables  and  stress  and 
velocity  gradients  are  the  secondary  variables.  This  unavoidably  complicates  understanding  and 
specification  of  boundary  conditions  for  these  variables,  and  it  can  be  counter  intuitive  when 
compared  to  solving  standard  stress  equilibrium  problems  which  do  not  involve  porous  media 
flow. 
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1.3.3  Interface  Modeling  Complexity.  The  high  order  element  and  mixed  element 
formulation  further  complicate  interface  modeling.  This  is  a  result  of  well-known  anomalous 
behavior  of  nodal  point  reaction  forces  in  high  order  elements.  A  simple  three-dimensional 
example  illustrates  this  point  in  Figure  1-3.  Though  the  computed  reaction  forces  are  in 
equilibrium  with  the  applied  load,  the  comer  node  reactions  are  in  fact  in  the  wrong  direction. 
This  anomalous  behavior  can  confuse  an  interface  computational  algorithm  that  imposes 
contact/release  conditions  on  node  points  along  the  plow/soil  interface  during  problem  solution. 

1.3.4  Accuracy  of  Constitutive  Models.  In  a  preceding  project,  advanced  finite  element 
technology  was  applied  to  a  geotechnical  problem  involving  penetration  of  rock  anchors  to  help 
resolve  discrepancies  between  field  observation  and  analytical  prediction  (Hoge  and  Shugar, 
1992).  DYNA3D  was  employed  in  this  study;  and  is  an  explicit  finite  element  program 
particularly  suitable  for  highly  transient  dynamic  problems  (Whirley  and  Hallquist,  1991).  This 
investigation  concluded  that  modeling  nonlinear  material  behavior  of  rock  was  the  primary 
challenge  and  that  advanced  constitutive  models  were  required  for  rock.  It  is  not  unexpected  that 
correctly  modeling  the  nonlinear  behavior  of  saturated  seafloor  soil  is  of  key  importance  in  the 
present  study  of  plow/soil  interaction.  Thus,  a  requirement  for  more  advanced  soil  models  may 
well  be  anticipated  in  the  long  term  for  solution  of  this  problem.  The  following  discussion 
provides  a  short  status  report  on  what  has  been  done  at  NFESC  in  this  regard. 

The  Office  of  Naval  Research  (ONR)  sponsored  the  development  of  modem,  third- 
generation  soil  constitutive  laws  for  use  in  numerical  modeling  of  geotechnical  systems  which 
was  conducted  at  NFESC,  the  University  of  California  at  Davis,  and  at  Stanford  University,  all 
under  the  auspices  of  the  NFESC’s  6.1  project,  Structural  Modeling.  Among  the  achievements  in 
computational  mechanics,  this  research  resulted  in  a  bounding  surface  plasticity  constitutive 
model  for  clay  soils  which  was  validated  in  a  fairly  wide  array  of  experiments  (Dafalias  and 
Herrmann,  1986;  Herrmann,  1987).  A  fundamental  experimental  study  on  the  behavior  of  sand 
was  also  completed  (Kutter,  Chen,  and  Shen,  1994).  In  this  study,  a  large  number  of 
experimental  validation  tests  on  a  medium  to  fine  grain  sand  were  conducted  to  provide  an 
experimental  data  base  for  development  and  validation  of  a  proposed  bounding  surface 
hypoplasticity  model  for  sand  (Wang,  Dafalias,  and  Shen,  1990)  and  third-generation  soil  model 
development  efforts  conducted  elsewhere.  Initial  attempts  to  validate  the  model  with  earthquake 
field  test  data  were  successful  (Shen  and  Li,  1991).  However,  initial  validation  results  were 
limited  to  one-dimensional  soil-column  models,  and  measured  pore  pressure  data  were  below 
critical  liquefaction  levels  and  therefore  insufficiently  challenged  the  model.  Further  validation 
efforts  involving  laboratory  data  including  geotechnical  centrifuge  experiments,  were  hampered 
by  numerical  inefficiency  and  also  revealed  shortcomings  in  the  physical  behavior  of  the 
proposed  hypoplasticity  model.  For  example,  the  hypoplasticity  sand  model  assumed  that  the  soil 
failure  envelop  for  undrained  conditions  is  the  same  as  that  for  drained  conditions,  but  laboratory 
data  indicated  otherwise.  The  experiments  further  showed  that  the  shape  of  the  phase 
transformation  surface  is  different  from  the  shape  of  the  failure  surface  and  suggested  that  the 
model  should  accommodate  this  difference  as  well. 

These  theoretical  and  numerical  problems  combine  to  place  at  risk  successful  validation 
of  the  hypoplasticity  model  for  sand  as  currently  proposed.  Additional  research  on  constitutive 
models  for  granular  materials  was  recommended.  The  numerical  problem  is  currently  being 
addressed  at  U.  C.  Davis  (Holland,  1996)  and  progress  is  being  made  by  analyzing  different 
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numerical  integration  algorithms  for  solving  the  so-called  incremental  rate  equations  (for 
evaluating  effective  stress  at  the  Gauss  points).  These  algorithms  are  central  to  efficient 
implementation  of  constitutive  models  in  general.  Residual  basic  research  exists  on  these 
problems  at  this  time,  and  no  additional  research  is  currently  programmed  to  address  either 
problem  under  ONR’s  basic  research  program.  In  spite  of  this,  a  recent  National  Science 
Foundation  (NSF)  sponsored  study  (Arulanandan  and  Scott,  1993)  is  encouraging,  and  has 
affirmed  the  potential  of  advanced  constitutive  models  and  their  numerical  implementation  (for 
liquefaction  problems)  in  geotechnical  engineering. 
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ANOMALOUS  BEHAVIOR  OF  HIGHER  ORDER  ELEMENTS 
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Figure  1-3.  Anomalous  behavior  of  higher  order  elements. 


2.0  ALTERNATIVE  MODELS  OF  STEADY-STATE  PLOW/SOIL  INTERACTION 


An  analysis  of  the  physical  processes  attending  plow/soil  interaction  was  an  early 
imperative  in  support  of  the  development  of  a  numerical  model.  An  engineering  services  contract 
was  prepared  for  U.C.  Davis.  The  statement  of  work  is  included  as  Appendix  B  and  may  be 
consulted  for  details. 

In  the  sequel,  two  specialized  approaches  to  plow/soil  interaction  modeling  were 
developed  at  U.C.  Davis  which  both  complement  and  support  the  present  modeling  task 
employing  commercial,  general-purpose  finite  element  technology.  The  contract  was  completed 
in  the  fall  of  1994  and  resulted  in  two  published  reports: 

1.  Kutter,  B.  L.  and  T.  Voss  (1995).  “Analysis  of  data  on  plow  resistance  in  dense, 
saturated  cohesionless  soil,”  CR  95.004-SHR,  Naval  Facilities  Engineering  Service 
Center,  Port  Hueneme,  CA  93043-4328,  February  1995. 

2.  Herrmann,  L.  R.  and  J.  Mello  (1994).  “Investigation  of  an  alternative  finite  element 
procedure:  a  one-step,  steady-state  analysis,”  CR  95.001-SHR,  Naval  Facilities 
Engineering  Service  Center,  Port  Hueneme,  CA  93043-4328,  December  1994. 

These  reports  have  been  reviewed  and  a  summary  of  their  results  and  the  plow/soil  interaction 
models  which  they  describe  follows.  The  respective  reports  may  be  consulted  for  detailed 
accounts. 

2.1  Approximate  Analytical  Model  -  A  Semi-Empirical  Approach  (Kutter  and  Voss,  1995) 

The  dilatancy  of  dense  sand  can  cause  large  negative  pore  pressures  to  develop  during 
rapid  plow  penetration.  Experimentally,  this  results  in  a  nearly  linear  increase  in  the  plow 
resistance  with  plow  velocity  up  to  a  critical  velocity  at  which  cavitation  of  the  pore  fluid  takes 
place.  Dimensional  analysis  is  useful  for  description  of  the  relationship  between  plow  resistance 
and  velocity.  The  number  of  required  parameters  is  reduced  by  the  use  of  dimensionless 
variables.  Critical  state  soil  mechanics,  bearing  capacity  theory,  and  Darcy’s  Law  are  used  for 
guidance  in  the  selection  of  alternative,  physically-admissible,  dimensionless  groups  of 
parameters.  The  experimental  data  base  (Cable,  et  al.,  1993;  Girard  and  Taylor,  1994)  provided  a 
means  to  empirically  evaluate  an  otherwise  unknown  function  of  dimensionless  parameters  that 
arose  in  the  derivation  of  the  admissible  parameter  groups.  The  result  was  a  semi-empirical 
model  which  appears  to  correlate  with  experimental  data,  and  which  adequately  describes  the 
rate  of  change  of  plow  resistance  with  velocity.  In  the  empirical  approach,  all  of  the  effects  of 
assumptions,  omitted  parameters,  etc.,  whether  or  not  significant,  are  assumed  to  be  lumped  into 
the  model  error,  —  the  distance  between  the  prediction  and  the  data.  Figure  2.1  schematically 
illustrates  this  error  along  with  semi-empirical  model  parameters. 

The  semi-empirical  plow/soil  interaction  model  is  based  on  existing  bearing  capacity 
equations  used  in  geotechnical  engineering  of  soil  foundations  (for  example,  see  Wu,  1967). 
These  equations  are  fundamentally  based  on  Prandtl  theory  of  slip  lines  formed  in  a  rigid  plastic 
material  where  bearing  stresses  are  maximum  and,  which  in  the  plow/soil  model,  emanate 
outward  into  the  soil  from  the  nose  of  the  plow.  The  model  does  account  for  three-dimensional 
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flow  and  deformation  of  soil  about  the  blade,  as  well  as  rake  angle  of  the  blade.  However,  the 
model  omits  parameters  not  included  in  the  data  base.  For  example,  the  model  assumes  that  side 
friction  on  the  blade  is  negligible  compared  to  bearing  resistance  at  the  nose  of  the  blade,  and  it 
has  no  parameter  for  blade  tip  included  angle.  Instead,  it  assumes  a  blunt  tip  thereby  ignoring  tip 
friction  as  well.  Friction  and  blade  sharpness  are  accounted  for  in  the  finite  element  numerical 
model  approach. 

The  important  properties  of  the  soil  which  appear  to  govern  the  mechanics  of  plow/soil 
interaction  are  listed  as  follows,  with  friction  angle  and  permeability  being  the  two  most 
important: 

1 .  Friction  angle,  0 

2  Initial  void  ratio,  eQ 

3.  Critical  state  void  ratio,  ecrjt 

4.  Permeability,  k 

5.  Vapor  pressure  of  pore  fluid,  Ucav 

6.  Ambient  pressure  of  pore  fluid,  uq 

The  important  geometric  and  kinematic  parameters  of  the  thin  blade  plow  which  govern 
its  resistance  to  steady-state  motion  according  to  the  present  semi-empirical  model  were  defined 
as  follows  (also  refer  to  Figure  2-1). 

1.  Length,  L 

2.  Thickness,  B 

3.  Depth,  H 

4.  Rake  angle,  0 

5.  Velocity  of  penetration,  v 

A  shear  zone  with  its  volume  assumed  to  be  proportional  to  the  cross-sectional  area  of  the 
blade  is  hypothesized.  It  comprises  a  region  of  uniform  volumetric  strain,  sv  ave,  which 
approximates  the  volume  of  actual  soil  otherwise  undergoing  very  nonuniform  dilatation  strain. 
This  volume  is  determined  empirically,  through  comparisons  of  the  model  with  experimental 
data  It  is  also  a  function  of  initial  and  critical  void  ratios  as  discussed  in  the  aforementioned 
(Kutter  and  Voss,  1995)  report.  Its  purpose  is  to  provide  for  the  linear  increase  in  plow  resistance 
with  velocity  observed  in  the  data. 

A  summary  of  the  procedure  for  estimating  plow  resistance  according  to  the  proposed 
semi-empirical  model  follows: 

1 .  Determine  cavitation  or  vapor  pressure  of  pore  fluid,  ucav 

2.  Determine  the  maximum  plow  resistance  by  the  smaller  of  either 
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Resistance 


Figure  2-1.  Definition  of  variables  for  semi-empirical  model  of  plow/soil  interaction 
(Kutter  and  Voss,  1995). 
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or. 


Q  max  -  (u0  ■  ucav)  NqBH 
Q  max  =  cu  NcBH. 


where  Nq  (a  function  of  friction  angle  -  0)  and  Nc  are  conventional  bearing  capacity 
factors,  respectively,  for  surcharge  and  for  cohesion,  and  Cu  is  the  undrained  shear 
strength  (Wu,  1967). 

3.  Estimate,  or  assume  negligibly  small,  plow  resistance  at  zero  velocity,  Q0 

4.  Estimate  rate  of  change  of  plow  resistance  with  velocity,  AQ/Av,  using 

AQ/Av  =  (y-w  NqBH  /k)  *  fQ(sy  ave>  B/H,  L/H,  q) 
where  yw  is  the  unit  weight  of  water. 

5.  Estimate  the  critical  velocity,  i.e.,  velocity  corresponding  to  maximum  plow 
resistance  using 


vcrit —  (Qmax  ■  QqV AQ/Av 


Values  for  the  unknown  soil  function  fQ  were  determined  empirically  using  the 
experimental  data  base.  They  range  between  0.0048  and  0.0122  for  soil  #1  (medium-dense  sand) 
and  between  0.0047  and  0.0071  for  soil  #2  (dense  sand),  depending  on  the  geometry  of  the  blade. 
This  function  also  depends  on  three-dimensional  aspects  of  soil  deformation,  and  it  may  also  be 
possible  to  evaluate  it  from  the  results  of  three-dimensional  numerical  models,  as  well  as  from 
experimental  data. 

It  was  concluded  that  the  important  soil  mechanics  features  of  plow/soil  interaction  are: 

1 .  Resistance  is  predominately  due  to  bearing  on  the  leading  edge  of  the  plow. 

2.  Shear  distortion  of  dense  sand  around  the  plow  causes  dilation,  and  the  rate  of  dilation 
is  a  function  of  the  rate  of  distortion,  which  is  proportional  to  the  velocity  of  plow 
penetration. 

3.  The  tendency  to  dilate  causes  negative  water  pressures  (relative  to  ambient  water 
pressure)  to  develop  in  the  shear  zone. 

4.  Negative  pressures  draw  water  into  the  shear  zone  at  a  rate  which  is  proportional  to 
the  change  in  pore  pressure  and  the  permeability  of  the  soil.  The  flow  of  water  into 
the  shear  zone  permits  dilation  of  particles  to  occur  and  relieves  the  negative  pore 
pressure. 
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5.  The  negative  pore  pressures  result  in  increased  effective  stresses  in  the  region  of  the 

deformation  of  the  soil. 

6.  Increased  effective  stress  results  in  increased  frictional  resistance  to  plow  penetration. 

Therefore,  plow  resistance  increases  as  velocity  of  penetration  increases. 

7.  The  maximum  magnitude  of  the  negative  pore  pressure  is  limited  by  either: 

a.  Cavitation  of  the  pore  fluid  when  the  absolute  pore  pressure  reaches  the  vapor 
pressure,  or 

b.  Achievement  of  a  critical  state  condition  in  the  shear  zone.  For  a  given  value  of 
void  ratio,  there  is  a  critical  value  of  the  effective  stress  at  which  the  soil  may 
shear  without  further  change  in  effective  stress,  pore  pressure,  or  shear  stress. 

Recommendations  were  also  given  for  additional  material  property  testing  and  for 
additional  laboratory  plow/soil  interaction  testing. 

2.2  Steady-State  Finite  Element  Analysis  Model  -  A  Rational  Model  Approach  (Herrmann 
and  Mello,  1994) 

A  steady-state  finite  element  analysis  procedure  that  determines  the  stresses  in  the  soil 
and  the  required  driving  force  to  pass  a  blade  through  soil  at  constant  velocity  has  been 
developed.  This  type  of  analysis  is  shown  to  be  a  potentially  attractive  alternative  to  a  multistep 
transient  analysis  (using  a  commercial  finite  element  program)  when  the  latter  becomes 
computationally  expensive.  A  steady-state  analysis  allows  time-dependence  to  be  expressed  as 
space-dependence,  and  thus,  a  solution  occurs  in  a  single  step.  The  feasibility  of  this  type  of 
analysis  was  demonstrated.  It  was  determined  that  more  investigation  into  soil  failure  conditions 
at  the  plow  tip  and  an  appropriate  soil  plasticity  model  for  those  conditions  are  required  before  a 
production  three-dimensional  single-step  (or  a  multistep)  computer  program  is  developed. 

A  plane  strain  formulation  for  steady-state  analysis  of  the  plow  soil/interaction  problem  is 
illustrated  in  Figures  2-2a  and  2-2b.  The  region  of  soil  which  is  depicted  in  plan  is  a  horizontal 
plane  which  has  unit  thickness  in  the  vertical  direction.  This  formulation  clearly  ignores  the 
proximity  of  the  free  surface.  The  dimensions  in  inches,  and  the  size  of  the  finite  element  model 
are  indicated  in  Figure  2-2c,  for  a  typical  analysis  with  a  steady-state  velocity  of  five  feet  per 
second. 

Results  from  the  model  predicted  plow  resistance  to  be  a  linear  function  of  plow  velocity 
as  has  been  observed  experimentally.  Further,  a  rather  dramatic  dependence  on  plow  shape,  i.e., 
tip  sharpness,  and  the  apparent  importance  of  including  inertia  forces  were  demonstrated  as 
shown  in  Figure  2-2d.  It  is  noted  that  soil  inertia  forces  are  explicitly  neglected  in  the  previous 
semi-empirical  model  for  it  has  no  mass  density  parameter.  It  was  assumed  that  they  were 
negligible. 
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(a)  Problem  configuration. 


(b)  Configuration  for  analysis. 


Figure  2-2.  Preliminary  single-step  steady-state  finite  element  model  of  plow/soil 
interaction  (Herrmann  and  Mello,  1994). 
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(c)  Deformed  mesh  for  plow  problem,  v0  =  5  fps. 


In  the  deformed  mesh  shown,  the  blade  is  represented  by  a  prescribed  rigid  surface  which 
is  “moving  with  steady-state  velocity”  along  the  centerline  of  the  solution  region,  otherwise  not 
visible  (due  to  the  expanded  scale  of  the  figure).  A  typical  result  for  soil  shear  stress  is  shown  in 
Figure  2-2 e,  where  it  is  clear  that  the  end  boundaries  of  the  model  are  sufficiently  distant  to  have 
not  affected  the  solution  region  near  the  plow.  This  isolation  is  more  difficult  to  achieve  with 
multistep  transient  numerical  models  because  it  results  in  greater  computational  effort. 

A  nearly  incompressible,  linear  viscoelastic  material  model  for  soil  was  used  in  this 
demonstration.  No  soil  inelasticity  or  pore  water  flow  was  considered  in  the  formulation  of  the 
model  and  stresses  are  total  stresses  not  effective  stresses.  However,  substantial  development  of 
appropriate  boundary  conditions  (including  derivative  boundary  conditions  and  plow-soil 
interface  modeling  conditions,  including  friction)  were  completed  to  accommodate  simulation  of 
conditions  imposed  in  the  experimental  test  apparatus  (Figure  1-1  a). 

Of  primary  significance  was  the  inexpensive  nature  of  the  steady-state  finite  element 
analysis.  The  model  consisted  of  2,800  degrees  of  freedom  and  required  only  a  matter  of  minutes 
to  solve  on  an  Intel  486-based  personal  computer.  Only  5  to  10  iterations  were  required  for 
convergence,  and  no  time-stepping  is  required. 

From  this  preliminary  study  several  conclusions  were  drawn.  Because  of  the  very  few 
iterations  (with  no  time-stepping)  required,  it  should  be  entirely  feasible  to  perform  a  three- 
dimensional,  steady-state,  inelasticity  study  of  the  plow/soil  interaction  problem  including 
parameter  studies  for  different  plow  shapes  and  plowing  velocities  (implementation  of  inelastic 
behavior  is  ongoing  at  U.  C.  Davis).  Thus,  the  single-step,  steady-state  analysis  procedure  offers 
a  viable  alternative  to  a  multistep,  transient  analysis  of  the  plow/soil  interaction  problem. 
However,  because  no  suitable  three-dimensional  steady-state  finite  element  program  currently 
exists  (so  far  as  is  known),  it  would  seem  that  the  best  way  toward  a  numerical  solution  is  to  use 
an  available  commercial  program,  to  perform  a  multistep,  transient  analysis. 

Two  possible  approaches  to  a  three-dimensional  nonlinear  plow/soil  interaction  model 
based  on  a  steady-state  finite  element  algorithm  were  suggested  along  with  a  few  tasks  necessary 
to  eliminate  obvious  deficiencies  in  the  current  software  (also,  see  Appendix  A  for  a  related 
discussion  involving  contingency  modeling  plans). 

1 .  The  present  two-dimensional  code  at  U.  C.  Davis  which  has  the  single-step,  steady- 

state  algorithm  can  be  extended  by  adding  the  following  capabilities: 

a.  Three-dimensional  finite  elements 

b.  Effective  stress  and  pore  water  flow  formulation 

c.  Nonlinear  plasticity  soil  material  models 

2.  NFESC’s  three-dimensional  nonlinear  finite  element  code  Geneseos  could  be 

modified  to  include  the  following  capabilities: 

a.  Data  structure  change  to  accommodate  the  steady-state  algorithm 

b.  Merge  existing  effective  stress  soil  model  with  new  data  structure 
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c.  Interface  element  and  slide  lines 


d.  Large  deformation 

In  either  case,  substantial  software  development  is  indicated.  However,  its  cost  can  be  reduced  if 
it  is  designed  and  written  in  conjunction  with  universities. 

Before  any  attempt  is  made  to  produce  a  steady-state  finite  element  production  program 
and  analysis  procedure,  several  items  require  further  investigation.  The  failure  condition  in  the 
soil  at  the  tip  of  the  plow,  i.e.,  the  state  of  the  soil  as  it  is  cut  by  the  plow  tip  to  form  a  new 
surface,  is  not  well  understood  and  is  an  area  that  should  receive  further  theoretical  and 
experimental  study.  (This  information  is  also  required  for  a  rigorous  multistep,  transient  analysis 
of  the  plow  problem.)  The  effect  of  using  a  realistic  plasticity  model  for  the  soil  on  the  rate  of 
convergence  of  the  iteration  process  in  the  single-step  model  must  be  studied.  Finally,  porous 
media  flow  must  be  coupled  into  the  analysis. 

Of  these  recommendations,  only  the  computational  plasticity  issue  is  currently  being 
addressed.  Moreover,  in  this  context,  a  finite  element  formulation  for  handling  strain  softening  in 
the  presence  of  inelastic  behavior  has  been  completed  under  Navy  6.1  funding  support  (Simo, 
1994;  Simo,  Oliver  and  Armero,  1993).  This  work  resulted  in  a  new  approach  to  the  analysis  and 
simulation  of  localization  arising  in  inelastic  solids  that  exhibit  strain-softening  response,  and  it 
has  been  extended  to  the  finite  deformation  regime.  The  technique  leads  to  the  systematic 
construction  of  numerical  models  that  completely  eliminate  the  strong  mesh  dependence 
exhibited  by  conventional  treatments  of  such  problems.  It  has  been  recently  extended  to  two- 
dimensional  numerical  simulations  (Armero,  1995),  and  warrants  further  research  in  support  of 
numerical  models  in  soil  mechanics  and  geotechnical  engineering  in  general. 
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(e)  Contour  plot  of  shear  stress  (txy)  in  soil,  v0=5  fps 


3.0  MULTISTEP  TRANSIENT  NUMERICAL  MODEL  OF  STEADY-STATE 
PLOW/SOIL  INTERACTION 

A  description  of  the  multistep  plow/soil  interaction  numerical  model  developed  in  this 
study  is  provided  in  the  following  paragraphs.  Since  the  problem  is  characterized  by  steady-state 
motion,  the  distinction  is  drawn  between  the  previous  single-step  model  approach  proposed  by 
Herrmann  and  Mello  (1994)  and  the  conventional  multistep  approach  employed  in  the  following 
development. 

3.1  Finite  Element  Formulation  for  Porous  Media 

ABAQUS  requires  that  a  high  order  finite  element  be  used  for  analysis  of  porous  media. 
These  elements  are  referred  to  as  CPE8P  and  C3D20P,  respectively,  for  two-  and  three- 
dimensional  applications.  The  two-dimensional  element  is  an  8-node  bi-quadratic  displacement 
and  bi-linear  pore  pressure  element,  and  the  three-dimensional  element  is  a  20-node  tri-quadratic 
displacement,  tri-linear  pore  pressure  element.  This  compares  with  4  and  8  nodes,  respectively, 
for  low  order  finite  element  formulations.  Use  of  the  high  order  elements  requires  managing 
larger  data  files  due  to  the  existence  of  midside  nodes  in  addition  to  comer  nodes.  Although,  high 
order  elements  are  generally  regarded  as  more  efficient  where  the  solution  involves  large 
gradients,  in  the  context  of  large-scale  computation,  low  order  elements  are  generally  favored  for 
modeling  simplicity. 

The  finite  element  formulation  for  porous  media  is  a  mixed  formulation;  displacements 
and  pore  pressures  are  primary  variables,  and  stress  and  pore  water  velocity  gradients  are 
secondary  (or  conjugate)  variables.  Corresponding  specification  of  boundary  conditions  for  these 
primary  and  secondary  variables  can,  at  first,  be  counter  intuitive.  For  example,  on  a  symmetry 
plane,  typically  the  displacement  normal  to  the  plane  is  constrained  to  zero.  However,  it  would 
be  incorrect  to  also  constrain  the  pressure  in  this  way,  because  this  simultaneously  implies  that 
flow  occurs  through  the  plane  in  obvious  violation  of  the  flow  boundary  condition.  Instead,  the 
pressure  variable  is  treated  as  a  free  variable  to  coerce  a  zero  flow  condition  normal  to  a 
symmetrical  boundary. 

Moreover,  higher  order  element  formulations  complicate  specification  of  interface 
models  involving  contact  forces.  This  is  a  result  of  well-known  anomalous  behavior  of  nodal 
point  reaction  forces  in  high  order  elements  (Figure  1-3).  Low  order  elements  do  not  have  this 
complication  and  tend  to  work  more  easily  with  interfaces. 

3.2  Basic  Plow/Soil  Interaction  Model  Concepts 

The  design  of  a  three-dimensional  plow/soil  interaction  model  has  been  considered.  The 
basic  model  concept  is  highly  constrained  in  that  there  are  no  known  alternative  approaches 
possible  using  ABAQUS/Standard.  In  discussions  with  technical  support  personnel  at  HKS,  Inc., 
it  was  discovered  that  the  program  does  not  provide  directly  for  contact  problems  in  porous 
media  that  involve  large  amounts  of  sliding  at  an  interface.  However,  its  overall  modeling 
capability  is  sufficiently  robust  that  the  following  proposed  approach  is  feasible  in  the  case  of 
either  two-  or  three-dimensional  modeling  of  the  plow/soil  interaction  problem. 
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A  sketch  of  the  concept  for  the  three-dimensional  case  is  presented  in  Figure  3-1.  The 
blade  is  modeled  with  rigid  body  elements.  The  soil  is  required  to  be  modeled  with  high  order, 
20-node  elements  having  the  coupled  equilibrium/porous  media  flow  formulation  (C3D20P). 
Between  the  blade  and  soil,  two-dimensional  interface  elements  are  required  to  model  the 
penetration  process.  However,  the  20-node  elements  have  the  aforementioned  anomalous 
behavior  involving  comer  node  forces.  Comer  node  reaction  forces  are  invariably  incorrect  in 
both  magnitude  and  direction  and  these  are  the  crucial  contact  forces  that  determine  the 
resistance  of  the  plow.  Their  usage  in  conjunction  with  contact  algorithms  can  potentially  cause 
havoc  with  solution  convergence,  and  thus  these  elements  are  generally  not  recommended  in 
contact  problems  (though  an  interface  element  is  available  for  use  with  them  in  ABAQUS). 
However,  by  also  specifying  the  multipoint  constraint  (MPC)  feature  of  ABAQUS  at  the 
interface,  the  high  order  soil  elements  interpolate  displacements  linearly  on  the  interface,  as  is  if 
they  were  low  order  elements  to  begin  with.  The  result  is  that  only  comer  nodes  are  active,  while 
midside  nodes  are  suppressed,  along  with  the  interface. 

'  By  activating  the  MPC  feature,  comer  node  reaction  forces  can  be  expected  to  behave 
more  smoothly.  Four-node,  vice  eight-node,  interface  elements  (IRS4)  are  then  placed  between 
the  rigid  plow  elements  and  the  deformable  soil  elements  to  maintain  compatibility  with  the 
constrained  20-node  soil  elements.  Subsequent  specification  of  frictional  behavior  always 
complicates  interface  models,  but  it  should  not  necessarily  be  more  complicated  because  of  this 
approach. 

Development  and  validation  of  a  three-dimensional  plow/soil  interaction  model  entails 
the  generation  of  a  finite  element  mesh,  presumably  in  which  the  basic  geometry  and  boundary 
conditions  of  the  mesh  would  be  selected  to  replicate  the  laboratory  test  apparatus  (Figure  1-la) 
as  closely  as  possible.  (ABAQUS  PRE,  an  automatic  pre-processor  with  a  graphics  user 
interface,  is  now  available  with  Version  5.5  to  assist  in  this  future  task).  For  the  most  part,  a 
coarse  mesh  is  anticipated.  However,  along  the  plow/soil  interface  a  finer  mesh  is  required  to 
capture  strong  deformation  gradients.  This  would  be  a  critical  consideration,  as  formidable 
compute  times  and  data  file  sizes  would  have  to  be  balanced  with  finite  element  mesh  size  and 
other  modeling  considerations  to  manage  a  three-dimensional,  multistep  transient  analysis  study. 

3.3  Basic  Two-Dimensional  Model 

As  discussed  for  the  single-step  steady-state  analysis  demonstration,  the  multistep 
transient  analysis  also  considered  a  plane  strain  formulation  in  the  horizontal  plane  just  beneath 
the  surface  of  the  soil.  The  assumption  of  a  plane  strain  state  of  stress  does  not  hold  up  well 
because  the  soil  stress  state  is  understood  to  be  three-dimensional  according  to  analysis  of 
experimental  data  from  laboratory  studies  (Cable  et  al.,  1993,  and  Girard  and  Taylor,  1994). 
These  data  show  that  the  stress  distribution  through  the  depth  is  not  uniform.  Further,  this  is  not 
surprising  given  the  proximity  of  a  ffee-surface.  Likewise,  a  plane  stress  assumption  is  also 
incorrect.  ABAQUS  supports  the  plane  strain  assumption  with  its  CPE8P  element,  which  is  the 
porous  media  counterpart  of  the  standard  high  order  quadrilateral  element  for  plane  strain,  the 
CPE8  element.  However,  there  is  no  corresponding  porous  media  element  counterpart  for  the 
standard  high  order  plane  stress  element,  CPS8.  Therefore  a  plane  strain  formulation,  as  depicted 
in  Figure  3-2,  was  pursued  to  demonstrate  the  numerical  model  approach  to  plow/soil  interaction. 


-  Infinity 


Figure  3-2.  Plane  strain  formulation  of  plow/soil  interaction. 
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The  actual  soil  resistance  force  Fr  has  been  observed  to  act  well  below  the  center  of  the 
test  article  (plow  blade)  in  the  laboratory  data,  so  it  follows  that  the  shape  of  the  actual  stress 
profile  would  be  roughly  linear  or  quadratic  (as  shown)  with  depth  from  the  free  surface.  The 
plane  strain  formulation  ignores  this  observed  profile  in  favor  of  a  uniform  stress  profile  (not 
shown)  extending  to  +  and  -  infinity  in  the  z-direction.  The  corresponding  plane  strain  solution 
domain  consists  of  a  unit  slice  of  this  system,  as  shown,  with  the  coordinate  system  and  a 
reaction  boundary  support  located  at  the  ends.  The  direct  strain  component  in  the  z-direction  is 
zero  in  the  solution  domain. 

This  two-dimensional  analysis  was  justified  on  a  preliminary  basis  due  to  the  anticipated 
numerical  modeling  complexity  of  the  plow/soil  interaction  problem.  For  example,  it  would  be 
necessary  to  model  soil  plasticity  behavior,  pore  water  flow,  large  deformations,  large  strain, 
steady-state  motion  of  rigid  surfaces  involving  substantial  penetration  of  a  rigid  blade,  and  a 
frictional  contact  interface  with  surrounding  soil.  Dealing  with  this  complexity  is  best 
accomplished  deliberately  with  a  series  of  two-dimensional  models  before  attempting  a  three- 
dimensional  numerical  model  of  plow/soil  interaction.  Once  these  issues  are  resolved  it  is 
presumed  that  an  extension  to  three-dimensions  would  be  at  least  conceptually  straightforward, 
considerations  of  large  compute  times  aside. 

Test  conditions  in  the  experimental  apparatus  (Figure  1-la)  guided  the  overall 
development  of  the  numerical  model  so  that  despite  the  crude  nature  of  a  plane  strain  formulation 
of  this  problem,  some  comparison  of  calculated  and  measured  response  would  be  qualitatively 
feasible. 

The  corresponding  plane  strain  finite  element  model  is  shown  in  plan  in  Figure  3-3.  One 
plane  of  symmetry  is  exploited  as  indicated  by  the  symmetry  line.  The  origin  of  the  global 
coordinate  system  is  placed  at  the  lower  left-hand  comer  of  the  soil  region,  with  the  x-y-z 
reference  axes  labeled  1-2-3,  respectively.  The  soil  region  shown  is  a  6-  by  12-inch  rectangular 
region,  and  the  blade  is  6.5  inches  long  with  a  1/4-inch  half-width,  simulating  a  1/2-inch  thick 
blade.  The  plane  strain  model  has  unit  thickness.  The  finite  element  mesh  is  finely  graded  along 
the  intended  path  of  travel  to  capture  expected  high  gradients  of  effective  stress,  strain,  pore 
pressure,  and  flow  response  in  the  region  of  the  blade  tip  during  penetration.  The  CPE8P  porous 
media  elements  are  high  order,  two-dimensional  quadrilateral  elements,  and  in  combination  with 
this  mesh,  should  provide  sufficient  resolution  of  the  computed  response  of  the  soil.  To  minimize 
computation  time,  prescribed  element  aspect  ratios  deteriorate  with  distance  from  the  centerline 
and  primary  solution  region.  The  motion  of  the  blade  is  steady-state,  allowing  the  blade  to 
penetrate  the  soil  to  a  distance  of  9.6  inches  where  the  simulation  is  terminated.  This  distance  is 
sufficient  for  steady-state  conditions  in  the  soil  response  to  evolve  while  providing  for  complete 
penetration  of  the  6.5-inch  long  blade.  It  is  probably  insufficient  to  effectively  isolate  the 
solution  domain  from  the  effects  of  artificial  boundaries. 

The  primary  nodal  variables  are  denoted  Ul,  U2,  and  p;  they  are,  respectively,  the  x-  and 
y-displacements  at  (both  comer  and  midside)  node  points,  and  the  pore  pressure  at  comer  node 
points.  The  boundary  conditions  at  the  left  edge  of  the  model  impose  symmetry  conditions  on 
these  variables.  The  Ul  component  must  remain  nonnegative  to  satisfy  this  condition  and  the  U2 
component  must  remain  free.  That  is,  soil  must  be  allowed  to  displace  in  all  directions  except 
through  the  symmetry  line.  Perhaps  not  so  obvious  is  that  pore  pressure,  p,  must  also  remain  free 
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(i.e.,  undetermined)  so  that  its  conjugate,  flow  velocity,  will  be  correspondingly  constrained  to 
zero  in  the  direction  normal  to  the  symmetry  line. 

All  ABAQUS  contact  elements  are  referred  to  loosely  as  “interface”  elements.  These 
elements  are  intended  for  rigid-to-deformable  body  contact  in  two-  and  three-dimensions.  The 
specific  rigid  surface  contact  element  used  in  the  plow/soil  interaction  model  is  the  IRS21.  This 
one-dimensional  element  has  two  nodes  at  either  end  and  two  Gauss  points  coincident  with  the 
nodes.  It  allows  large  relative  motion  to  occur  between  contacting  bodies,  and  the  rigid  surface 
element  output  data  are  contact  pressures  (a  positive  quantity  reported  at  the  nodes).  Rigid 
surface  contact  elements  also  have  a  geometry  specification  option  for  defining  the  shape  of  a 
rigid  body  which  corresponds  to  the  contact  element.  These  shapes  have  a  single  node  located  at 
the  center  of  the  body  which  is  known  as  a  rigid  body  reference  node.  Since  rigid  bodies  may 
also  have  motion,  the  motion  of  this  node  must  also  be  appropriately  constrained.  The  IRS21 
contact  elements  are  located  on  the  surface  of  the  deformable  body  which  is  to  come  into  contact 
with  the  specified  rigid  body. 

To  implement  the  inequality  constraint  U1>0  and  to  provide  for  large  sliding  motion  to 
occur  along  the  symmetry  line,  this  line  is  regarded  as  rigid  surface  along  which  the  plow,  also 
regarded  as  rigid,  slides.  A  series  of  IRS21  contact  elements  are  therefore  specified  along  the 
symmetry  line,  each  of  which  corresponds  to  a  rigid  body  reference  node  at  the  center  of  the 
specified  rigid  surface  (extending  from  y  =  -1  to  y  =  13  inches).  Since  this  rigid  surface  must 
remain  fixed,  its  reference  node  is  completely  restrained.  The  rigid  surface  has  the  effect  of  a 
wall  through  which  soil  elements  may  not  penetrate  as  they  are  deformed  in  advance  of  the 
approaching  plow. 

Once  the  plow  has  arrived  and  displaces  soil  away  from  the  wall,  large  relative  motions 
and  contact  conditions  must  also  be  provided  for  at  the  interface  between  soil  and  the  plow.  To 
this  end,  a  second  series  of  IRS21  contact  elements  are  specified  along  the  symmetry  line,  one 
for  each  of  the  75  soil  elements.  Unlike  the  first  series,  which  are  linked  to  the  rigid  wall  surface, 
this  series  is  linked  to  the  rigid  body  reference  node  at  the  center  of  a  prescribed  rigid  surface 
having  the  shape  of  the  half  plow  blade,  as  shown.  This  rigid  body  reference  node  is  completely 
restrained  except  for  translation  along  the  centerline.  In  this  direction,  the  plow’s  reference  node 
is  constrained  to  move  at  a  steady-state  velocity.  The  corresponding  reaction  forces  of  interest 
which  are  reported  at  this  node  are  the  drawbar,  or  tow  force,  and  the  lateral  force  (RF1  and  RF2, 
respectively)  necessary  to  coerce  the  prescribed  steady-state  velocity  for  the  plow  through  the 
soil. 

Specification  of  the  interface  behavior  for  the  plow/soil  interaction  model  was  involved. 
(Version  5.4  ABAQUS/Standard  is  reported  to  be  substantially  improved  in  this  regard).  This 
behavior,  however,  may  be  thought  of  simply  as  a  parting  of  the  two  prescribed  layers  of  IRS21 
contact  elements  which  are  initially  co-linear  along  the  symmetry  line,  and  which  are  later 
separated  by  the  advancing  plow.  One  layer  sticks  to  the  wall  and  the  other  is  pried  away  and 
becomes  sandwiched  between  the  plow  blade  and  the  soil  while  the  blade  passes  by. 

Two  alternative  failure  modes  for  the  soil  are  supported  by  these  interface  models  along 
the  centerline.  They  are  indicated  in  Figure  3-4,  and  are  termed  a  displacement  failure  mode  and 
a  tearing  failure  mode.  In  either  case  the  plow  must  slide  along  the  centerline  and  the  soil  must 
slide  over  the  blade.  The  difference  is  in  the  conditions  under  which  the  soil  separates  from  the 
centerline  ahead  of  the  plow  tip. 
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Figure  3-3.  Two-dimensional  plane  strain  finite  element  model  of  plow/soil  interaction. 
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Symmetry  line 


Displacement  mode  l  Tearing  mode 


Figure  3-4.  Alternative  failure  modes. 
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Along  the  upper  boundary  of  the  plow/soil  interaction  model,  the  U2  displacement 
component  is  constrained  to  zero  to  provide  a  reaction  to  the  plow  draw  bar  force,  and  simulate 
rather  coarsely  the  distant  rigid  end  of  the  experimental  apparatus,  which  is  a  long  narrow  box 
filled  with  soil  (Figure  1-1  a). 

ABAQUS  has  no  cavitation  pressure  cutoff  capability,  and  drainage  paths  in  the  z- 
direction  cannot  exist  in  the  two-dimensional  model  so  negative  pore  pressures  tend  to  large 
values  as  the  soil  dilates.  This  can  lead  to  incorrect  prediction  of  negative  pore  pressure  beyond 
the  cavitation  pressure  associated  with  the  test  apparatus  (about  -14.5  psi).  To  mitigate  this 
difficulty,  the  pressure  variable,  p,  is  constrained  to  zero  along  the  edge  line  (x  =  6  inches) 
allowing  flow  normal  to  this  boundary  as  if  the  side  wall  was  porous.  The  prescribed 
displacement  boundary  conditions  along  this  edge  are  designed  to  simulate  a  rigid  side  wall,  i.e., 
U1  is  zero  and  U2  is  free. 

Soil  displacement  at  the  lower  edge  of  the  model  (y  =  0)  is  unconstrained  to  allow  soil 
displaced  by  embedment  of  the  rigid  blade  to  escape  confinement.  Upon  constraining  this 
displacement,  it  was  found  that  penetration  of  the  blade  was  prevented  due  to  nearly 
incompressible  behavior  in  the  soil.  Also,  soil  in  the  experimental  test  apparatus  is  displaced 
vertically  by  the  blade  rather  easily  due  to  the  presence  of  the  free  surface,  and  this  unconstrained 
boundary  at  the  lower  edge  partially  compensates  for  the  model  having  no  such  free  surface. 

3.4  Analysis  with  Preliminary  Two-Dimensional  Numerical  Models 

The  aforementioned  features  of  the  plow/soil  interaction  model  were  not  all  implemented 
at  once.  Two  preliminary  models  were  developed  and  studied  from  which  these  features  evolved. 

3.4.1  Two-Dimensional  Elastic  Analysis  Model.  A  preliminary  two-dimensional 
elastic  analysis  initially  demonstrated  that  penetration  of  the  blade  into  the  soil  could  indeed  be 
simulated  in  a  satisfactory  manner.  A  multistep,  steady-state  transient  analysis  using  eight-node 
quadrilateral  elements  (CPE8;  or  dry  soil  elements)  for  the  soil  was  conducted  for  this 
demonstration.  Linear  elastic  properties  prescribed  for  the  soil  were  Young’s  modulus  equal  to 
6,000  psi  and  Poisson  ratio  equal  to  0.2. 

This  demonstration  incorporated  most  of  the  aforementioned  interface  and  rigid  surface 
modeling  techniques  necessary  to  simulate  large  sliding,  and  constituted  a  reasonable  beginning 
for  the  plow/soil  interaction  model.  A  one-dimensional  high  order  rigid  surface  contact  element 
(IRS22)  designed  to  work  directly  with  the  two-dimensional  high  order  element  was  used.  The 
multipoint  constraint  condition  was  not  specified  for  the  soil  nodes  along  the  centerline  in  this 
initial  demonstration.  Constant  velocity  (1  ft/s  and  8  ft/s)  penetration  was  prescribed  for  the 
blade,  and  it  was  successfully  advanced  along  the  sliding  surface  into  the  soil  under  quasi-static 
conditions. 

Run  times  for  this  preliminary  model  were  from  3-1/2  to  4-1/2  hours  on  a  Sparc  2 
Workstation.  This  may  seem  long  for  a  linear  elastic  characterization  of  the  soil,  but  interface 
elements  are  essentially  nonlinear,  and  therefore  the  analysis  is  nonlinear.  During  a  computer 
run,  many  iterations  of  the  contact  algorithm  were  required  to  resolve  inconsistencies  between 
both  gaps/over-closures  and  forces  at  the  interfaces.  Furthermore,  nonlinear  geometric  behavior 
was  specified  to  provide  for  large  displacements  in  the  soil  near  the  plow  blade. 
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A  dynamic  analysis  of  the  plow/soil  interaction  model  could  not  be  performed,  since  no 
inertia  is  included  in  the  ABAQUS  formulation  for  transient,  coupled  soil-structure  interaction 
analysis.  Inclusion  of  inertia  forces  in  the  soil  is  crucial  for  earthquake  problems,  for  example, 
but  was  generally  considered  unimportant  to  the  present  problem  which  was  otherwise  regarded 
as  a  low  velocity  penetration  problem.  The  kinetic  energy  transferred  from  the  blade  to  the  soil 
correspondingly  was  assumed  to  be  negligible.  Kutter  and  Voss  (1995)  seem  to  agree  with  this 
assumption,  at  least  for  smaller  velocities,  while  Herrmann  and  Mello  (1994)  indicate  that  soil 
inertia  forces  may  be  important.  Nonetheless,  quasi-static  procedures  were  employed  to  start  the 
plow  from  rest  and  carry  it  through  the  soil  in  this  initial  two-dimensional  elastic  plow/soil 
interaction  model. 

3.4.2  Two-Dimensional  Elastic-Plastic  Analysis  Model.  In  this  analysis,  the  added 
complexity  of  material  nonlinearity  was  introduced  into  the  plow/soil  interaction  model.  The  soil 
was  modeled  as  an  elastic-plastic  solid.  This  preliminary  model  was  also  prompted  by  solution 
convergence  difficulty  with  a  Drucker-Prager  model  specification  for  dry  soil,  i.e.,  a  total  stress 
version  of  the  plow/soil  interaction  model,  which  was  the  next  model  studied  after  the  elastic 
analysis  described  above.  (The  Drucker-Prager  model  was  eventually  made  to  work  (see  Section 
3.5.1)).  As  before,  the  elastic  properties  were  Young’s  modulus  equal  to  6,000  psi  and  Poisson 
ratio  equal  to  0.2.  The  yield  strength  was  set  to  30  psi,  and  strain  hardening  was  defined  by 
specifying  an  ultimate  strength  of  either  800  or  900  psi  at  various  yield  strains,  beyond  which  the 
material  model  behaves  perfectly  plastic.  The  yield  strength  of  30  psi  was  representative  of  the 
failure  stress  for  a  dense  sand.  Four  different  strain-hardening  curves  were  tested,  only  one  of 
which  led  to  a  converged  solution,  as  indicated  in  Figure  3-5.  Convergence  of  the  plasticity 
algorithm  was  often  prevented  by  too  much  deformation  in  the  soil  near  the  shoulder  of  the  plow 
tip.  Successful  run  times  typically  required  about  seven  hours  for  the  two-dimensional  elastic- 
plastic  plow/soil  interaction  model. 

3.5  Drucker-Prager  Material  Model  Calibration 

This  task  involved  studying  the  response  of  the  Drucker-Prager  plasticity  material  model 
under  various  conditions  with  a  simple  finite  element  model  of  a  triaxial  test  specimen.  These 
conditions  included  use  of  this  model  with  dry  soil  and  with  pore  water  under  both  drained  and 
undrained  conditions,  and  with  both  associative  and  nonassociative  flow  laws.  Elastic  behavior  is 
also  modeled  by  a  linear  elasticity  model  (*ELASTIC)  for  greater  efficiency  compared  to  the 
optional  porous  elasticity  model  (*POROUS  ELASTIC).  Although  the  latter  might  be  more 
appropriate,  most  of  the  strain  in  the  neighborhood  of  the  plow  was  plastic  strain.  These  options 
and  those  which  are  discussed  below  for  controlling  the  elastic-plastic  behavior  of  the  soil  model 
are  prescribed  under  the  *  MATERIAL  command  in  ABAQUS  (see  Appendix  C). 

The  material  model  study  of  a  triaxial  specimen  was  carried  out  in  three  dimensions.  The 
same  finite  element  formulation  used  in  the  two-dimensional  plow/soil  interaction  model  was 
also  used  in  this  three-dimensional  model  of  the  triaxial  specimen.  Success  was  based  on  how 
well  the  measured  response  of  triaxial  soil  samples  was  replicated  using  the  Drucker-Prager  soil 
material  model  in  conjunction  with  the  three-dimensional  model  of  triaxial  test  specimens. 
Reasonable  agreement  with  measured  response  constituted  calibration  of  the  Drucker-Prager 
model.  A  graph  of  triaxial  test  data  showing  deviator  stress  versus  axial  strain  was  used  to 
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partially  calibrate  the  Drucker-Prager  model  directly.  Consequently,  the  calibrated  model 
reproduces  this  data  accurately,  within  an  arbitrary  number  of  piecewise  linear  segments. 

3.5.1  Modified  Drucker-Prager  Material  Model.  A  Modified  Drucker-Prager  material 
model  was  used  to  simulate  soil  in  the  plow/soil  interaction  study.  The  variables  and  features  of 
this  model  are  indicated  in  Figure  3-6.  It  is  specifically  more  suited  to  analysis  involving 
monotonic  loading,  such  as  for  limit  load  analysis  of  a  soil  foundation.  It  assumes  that  a  purely 
elastic  response  regime  exists  within  the  yield  surface,  shown  in  Figure  3-6a  in  principal  stress 
space,  beyond  which  some  of  the  material  deformation  is  not  recoverable  and  can  be  idealized  as 
being  plastic.  Yield  behavior  depends  on  hydrostatic  pressure,  so  that  the  material  becomes 
stronger  with  hydrostatic  confinement  stress  as  shown.  However,  since  the  yield  surface  is  not 
closed,  this  model  cannot  predict  compaction  well  in  consolidation  type  problems.  To  do  this  the 
modified  Cap  model  in  ABAQUS  would  be  required.  The  Modified  Drucker-Prager  model  also 
provides  for  isotropic  strain  hardening  and  strain  softening  of  the  material.  The  input  data  is 
described  in  Table  3-1. 

Inelastic  behavior  is  generally  accompanied  by  plastic  volume  change  for  sand.  Either  an 
associative  or  a  nonassociative  plastic  strain  flow  rule  is  allowed,  according  to  whether  the  angle 

of  dilatancy,  ¥,  equals  or  is  less  than  the  friction  angle,  (3.  As  indicated  in  Figure  3-6b,  in  the 
plane  of  deviatoric  stress,  q,  and  mean  pressure  stress,  p,  (i.e.,  the  meridional  plane)  in  the  former 

case,  the  direction  of  an  increment  of  plastic  strain,  deP^,  is  taken  to  be  normal  to  the  yield 
surface  and  consequently  no  volume  change  takes  place  during  material  loading.  In  the  latter 
case,  the  loading  direction  has  a  component  parallel  to  the  yield  surface  which  indicates  the 
amount  of  plastic  volume  change.  For  sands  in  general,  the  experimental  and  analytical  evidence 

is  clear  that  ¥  is  less  than  (3.  Thus,  a  nonassociative  flow  rule  was  specified  exclusively 
throughout  this  material  study,  and  for  the  plow/soil  interaction  model. 

3.5.2  Calibration  of  Stress-Strain  Behavior.  The  calibration  method  used  for  drained 
conditions  is  illustrated  in  Figure  3-7.  A  physical  triaxial  soil  sample  has  a  cylindrical  shape.  The 
method  of  loading  and  the  flow  of  pore  water  during  consolidation  in  this  case  are  indicated  in 
Figure  3-7a.  The  stress  state  and  response  of  the  sample  is  generally  regarded  as  homogeneous 
until  localized  failure  of  the  sample  occurs.  That  is,  the  soil  state  is  assumed  the  same  throughout 
the  sample  during  any  instant  of  loading  up  to  localized  failure.  To  model  the  triaxial  sample,  a 
2  by  2  by  2  cube  of  soil  at  the  center  of  the  triaxial  specimen,  as  shown  in  Figure  3 -7b,  is 
considered  to  be  a  representative  sample  of  the  specimen  undergoing  uniform  behavior.  A  finite 
element  model  of  this  cube  is  employed  rather  than  of  the  entire  cylindrical  sample.  Exploiting 
three  planes  of  symmetry,  the  analysis  is  accomplished  with  a  single,  unit  cube,  finite  element  as 
shown  in  Figure  3-7c.  Thus,  the  material  model  study  is  carried  out  in  an  expedient  way  which  is 
reasonably  accurate  for  calibration,  and  which  is  also  representative  of  the  finite  element 
formulation  which  is  to  be  employed  in  the  plow/soil  interaction  model. 

The  model  is  executed  in  two  successive  load  phases.  The  first  is  a  geostatic  phase 
whereby  initial  effective  stresses,  pore  pressures,  and  confinement  stresses  are  prescribed  and 
allowed  to  gradually  seek  equilibrium.  The  second  phase  simulates  the  triaxial  test  and  prescribes 
a  uniform  displacement  to  the  top  surface,  compressing  the  model  to  20  percent  axial  strain,  as 
shown  in  Figure  3-8;  deformation  of  the  unit  cube  is  indicated  in  Figure  3-8a,  and  the 
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corresponding  pore  pressure  distribution  on  the  surface  is  indicated  in  Figure  3-8b.  As  indicated 
in  (Figure  3-7c),  node  1  at  the  origin  is  completely  restrained,  and  also  the  path  of  fluid  flow  is 
along  the  major  diagonal  from  node  209  to  node  1.  In  these  data  the  volume  increase  due  to 
imposed  deviatoric  stress  and  dilation,  and  the  negative  pore  pressure  contours  and  pore  pressure 
gradients  on  the  cube  surface  are  consistent  with  flow  towards  the  origin. 

Stress-strain  triaxial  test  data  for  soil  #1  from  Cable  et  al.  (1993),  are  plotted  in  Figure 
3-9.  Experimental  strain-hardening  behavior  is  indicated  in  the  graph  of  deviator  stress  versus 
axial  strain  in  Figure  3 -9a.  The  corresponding  model  result  is  shown  (solid  line),  where  it  is  seen 
that  the  measured  behavior  was  replicated  by  the  calibrated  model.  Both  strain  hardening  and 
softening  are  clearly  required  in  the  Drucker-Prager  model  if  it  is  to  replicate  the  measured 
triaxial  test  behavior.  By  contrast,  a  simplified  (elastic-perfectly  plastic  material)  version  of  the 
Drucker-Prager  model  can  capture  the  residual  stress  behavior,  but  not  the  hardening  and 
softening  behavior. 

The  accompanying  volume  change  in  the  measured  triaxial  data  is  only  calibrated  well 
for  strains  up  to  a  certain  level,  after  which  the  model  over  predicts  volume  change.  Volume 
change  is  simulated  by  specification  of  a  nonassociative  flow  model.  The  calibrated  value  for  the 

dilation  angle,  V  =  21.5  degrees,  was  obtained  through  a  trial  and  error  process.  The  model 
begins  to  over  predict  volume  expansion  of  the  triaxial  test  specimen  for  axial  strains  larger  than 
about  1 1  to  1 5  percent,  as  shown  in  Figure  3-9b. 

The  strength  of  the  material  is  prescribed  by  a  fixed  value  for  friction  angle,  which  was 
determined  from  experimental  data  to  be  p  =  46.2  degrees.  It  was  prescribed,  along  with  the 
dilation  angle,  using  the  *DRUCKER  PRAGER  option  which  describes  the  yield  surface  for  the 
elastic-plastic  model.  That  these  values  are  not  allowed  to  vary  with  material  state,  introduces 
error  and  accounts  for  the  discrepancy  in  modeling  volume  expansion. 

Likewise,  the  amount  of  isotropic  strain  hardening  and  softening  in  the  Modified 
Drucker-Prager  model  is  not  allowed  to  vary  with  the  material  state,  and  otherwise  artificially 
controlled  by  specifying  a  value  for  the  cohesion  parameter,  d  (Figure  3-6b).  The  input  data  for 
hardening  and  softening  behavior  is  listed  in  Table  3-1,  and  graphed  in  Figure  3-10.  These  data 
are  the  plastic  part  of  the  triaxial  data  (Figure  3-9a). 
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Figure  3-5.  Elastic-plastic  soil  models  for  2-D  plane  strain  plow/soil  interaction  model 
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Figure  3-6.  Modified  Drucker-Prager  plasticity  model 
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Figure  3.7.  Drained  triaxial  test  simulation  approach  for  calibration  of  Drucker-Prager 
material  model. 
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3-8.  Triaxial  test  numerical  model  results 
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(b)  Pore  pressure  distribution  on  surface  of  cube. 
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(b)  Calibrated  volumetric  behavior. 
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Figure  3-10.  Yield  stress  versus  plastic  axial  strain  in  compression. 


Table  3-1 .  Input  Table  for  Calibrating  Strain  Hardening/Softening 
Behavior  in  Drucker-Prager  Model 


Deviator  Stress 
(psi) 

Plastic  Component  of  Axial  Strain 
(in.) 

0.1 

0.0 

14.0 

0.005 

21.0 

0.010 

25.0 

0.015 

27.0 

0.021 

27.0 

0.040 

15.0 

0.140 

3.6  Two-Dimensional  Effective  Stress  Model  with  Special  Handling  Features 

The  plow/soil  interaction  model  is  based  on  a  nonlinear  quasi-static  coupled 
deformation/flow  analysis.  A  multistep,  steady-state  transient  analysis  solution  is  employed 
using  a  simplified  version  of  the  aforementioned  Modified  Drucker-Prager  effective  stress  soil 
model. 

A  description  of  four  potential  modeling  deficiencies  in  the  problem  of  plow/soil 
interaction  follows: 

1 .  Inertia  forces  are  neglected  in  the  formulation  for  dynamic  soil-structure  interaction. 
It  is  not  clear  that  these  forces  are  significant  in  the  present  problem,  but  there  is  no 
way  to  know  without  this  capability. 

2.  Cavitation  in  the  fluid  phase  cannot  be  simulated.  This  was  initially  thought  to  have 
been  an  insignificant  limitation  because  the  experimental  data  showed  evidence  of 
pore  pressure  values  near,  but  below,  cavitation.  However,  negative  (or  tensile)  pore 
pressures  greater  than  the  cavitation  pressure  were  predicted  by  the  model  and  there 
was  no  satisfactory  way  of  preventing  these  high  negative  pore  pressures. 

3.  Problem  solution  time  is  very  large.  This  makes  it  difficult  to  anticipate  investigation 
of  model  parameters,  although  some  parametric  studies  were  nonetheless  pursued. 
The  required  CPE8P  soil  element  is  a  high  order  element  and  a  major  reason  for  large 
solution  time.  The  two-dimensional  model  required  about  one  week  to  run,  and  a 
three-dimensional  model  would  be  prohibitive  on  the  Sparc  2  workstation.  A  three- 
dimensional  model  would  most  certainly  require  substantial  super  computer  access 
time.  In  a  typical  computer  run,  hundreds  of  gap  over-closures  and  openings  occur  at 
interfaces,  causing  the  system  to  cut  back  on  the  increment  size  and  to  try  again  until 
the  impenetrability  constraints  have  been  satisfied  within  a  prescribed  tolerance. 
Occasionally,  negative  eigen  values  are  computed  and  the  system  warns  of  numerical 
instability,  but  the  system  generally  does  recover. 
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4.  The  Modified  Drucker-Prager  soil  material  model  is  overly  sensitive  to  tensile  values 
of  mean  pressure  which  inevitably  develop  in  small  isolated  regions  near  the  tip  and 
wake  of  the  blade.  Consequent  numerical  problems  would  cause  the  simulation  to 
crash.  To  overcome  this,  the  calibrated  hardening/softening  behavior  in  the  material 
model  had  to  be  discarded  in  favor  of  a  simplified  elastic-perfectly  plastic  behavior 
(Figure  3 -9a).  This  allowed  a  closed  form  algorithm  extant  for  automatic  calculation 
of  a  material  stiffness  to  be  invoked  to  avoid  numerical  problems  and  allow  the 
simulation  to  continue.  As  stated,  this  simplified  model  does  replicate  residual  stress 
behavior,  and  further  strains  should  be  large  at  the  interface  of  the  blade,  where  the 
material  is  expected  to  be  in  a  failure  state.  Unfortunately,  this  could  not  be  shown  by 
calculation  since  the  plow/soil  interaction  model  could  not  be  made  to  work  unless 
the  Modified  Drucker-Prager  material  model  was  first  simplified  as  discussed. 

3.6.1  Simplified  Drucker-Prager  Soil  Model.  The  plow/soil  interaction  model  worked 
satisfactorily  with  an  elastic  soil  model.  Further,  the  Drucker-Prager  constitutive  model  was 
satisfactorily  calibrated  using  a  single  three-dimensional  finite  element  having  the  same  basic 
formulation  as  the  two-dimensional  finite  element  employed  in  the  plow/soil  interaction  model 
(Figure  3-3).  The  elastic  model  study  and  the  material  model  calibration  are  recommended 
preliminary  steps  in  the  development  of  a  complex  model  for  plow/soil  interaction.  There  was  a 
reasonable  expectation  that  the  basic  plow/soil  interaction  model  in  combination  with  the 
calibrated  Modified  Drucker-Prager  soil  model  would  then  also  work  smoothly  enough  from  the 
beginning.  Unfortunately,  this  was  not  true,  due  to  a  numerical  problem  associated  with  a 
singularity  in  the  yield  surface  of  the  Drucker-Prager  model  which  did  not  manifest  itself  during 
the  rather  extensive  calibration  study. 

The  source  of  the  singularity  can  be  seen  in  Figure  3-6a,  for  example,  where  at  the  vertex 
of  the  yield  surface  in  the  tension  region,  the  gradient  of  the  yield  surface  is  undefined.  When  the 
mean  stress  at  any  point  in  the  soil  becomes  tensile  (negative)  conditions  are  created  such  that  the 
rate  equations  in  the  formulation  of  the  Drucker-Prager  model  will  not  converge  to  a  solution.  A 
tensile  stress  state,  as  represented  by  a  stress  point  on  the  tensile  part  of  the  yield  surface,  will 
generally  tend  to  gravitate  towards  this  vertex.  Since  the  formulated  rate  equations  include  terms 
involving  the  gradient  of  the  surface,  which  at  the  vertex  is  undefined,  the  solution  for  the  stress 
state  in  the  soil  will  not  converge.  This  often  happens,  for  example,  in  the  analysis  of  soils  with 
large  angles  of  friction  (46.2  degrees  in  the  present  study)  and  causes  these  stress-point 
integration  algorithms  to  become  numerically  unstable. 

ABAQUS  signals  plasticity  algorithm  warnings  during  the  computation  in  these 
instances,  and  nonconvergence  was  experienced  in  some  cases  even  for  very  small  or  initial 

penetration  of  the  plow  tip,  and  even  for  very  small  time-step  sizes  of  order  10  .Very  isolated 
and  very  minor  tensile  values  were  invariably  noted  in  contour  plots  of  mean  stress  just  ahead  of 
the  plow  tip,  and  also  sometimes  in  the  wake  of  the  plow  blade.  They  were  otherwise  dominated 
by  compression  throughout  the  soil  region.  Nonetheless,  these  insignificant  tensile  stresses  were 
apparently  the  source  of  the  difficulty. 

Several  unsuccessful  attempts  were  made  to  eliminate  the  possibility  of  tensile  stress  by 
precompressing  the  soil  region  using  a  combination  of  prescribed  boundary  conditions  for  stress 
and  displacement.  However,  the  difficulty  could  not  be  eliminated,  and  it  was  concluded  that, 
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because  of  the  singularity  problem  associated  with  usage  of  the  Drucker-Prager  material  model, 
the  distributed  version  of  ABAQUS,  Version  5.3,  was  not  capable  of  executing  the  plow/soil 
interaction  model  without  substantially  simplifying  the  calibrated  Drucker-Prager  material 
model. 

A  development  version  of  ABAQUS,  Version  5.4,  (for  which  a  special  UMAT  routine  (a 
user-provided  FORTRAN  subroutine  for  a  material  model)  for  the  Drucker-Prager  model  has 
been  written  at  HKS  (West),  Inc.)  proved  successful  in  executing  the  plow/soil  interaction  model. 
The  UMAT  routine  included  an  artificial  stiffness  to  stabilized  numerical  difficulty. 
Unfortunately,  it  was  found  that  the  output  data  files  from  Version  5.4  could  not  be  post- 
processed  by  the  distributed  Version  5.3  of  POST,  the  ABAQUS  post-processor,  so  that  no  data 
from  the  successful  simulation  was  available.  Further,  this  unreleased  version  of  ABAQUS  was 
generally  unavailable  outside  HKS,  Inc.,  at  the  time. 

The  course  pursued  was  to  simplify  the  Drucker-Prager  model  so  that  the  troublesome 
stress-point  algorithm  would  be  avoided  in  the  calculation  of  stress  at  the  Gauss  integration 
points  within  each  element.  This  would  allow  Version  5.3  to  be  used  successfully,  but  at  the 
expense  of  the  Drucker-Prager  model  calibration  effort  which  replicated  the  hardening  and 
softening  behavior  evident  in  the  laboratory  triaxial  test  data  for  the  sand  material. 

A  software  fix  for  the  vertex  problem  in  the  Drucker-Prager  model  exists  in  Version  5.3, 
but  only  if  the  model  is  simplified  by  eliminating  the  strain  hardening  and  softening  parts  of  the 
stress-strain  behavior  in  favor  of  a  simpler  elastic-perfectly  plastic  model  of  stress-strain 
behavior  (Figure  3-9a).  With  this  simplification,  a  closed  form  analytical  formulation  is 
substituted  for  the  more  complicated  numerical  stress-point  algorithm  when  computing  stress  at 
Gauss  points.  Clearly,  this  action  is  driven  by  convenience  and  is  somewhat  arbitrary  but  it 
obviates  numerical  convergence  difficulties  with  the  material  model.  Further,  an  argument  can  be 
made  that  the  residual  stress  in  the  triaxial  behavior  of  the  sand  material,  i.e.,  the  stress 
corresponding  to  very  large  strain,  is  the  important  part  of  the  material  behavior  in  the  present 
study.  This  is  because  the  material  next  to  the  plow  is  obviously  subjected  to  very  large  strain, 
and  thus  levels  of  bearing  stress  on  the  plow  are  likely  to  be  representative  of  residual  stress 
levels. 

In  spite  of  this  simplification,  the  friction  angle  which  accounts  for  soil  strength  and  the 
angle  of  dilation  which  accounts  for  volume  expansion  of  the  Drucker-Prager  material,  both  of 
which  derive  from  experimental  observation,  remain  intact  and  valid  for  the  material  model  used 
in  all  subsequent  simulations  of  the  plow/soil  model.  Further,  though  the  solid  phase 
characterization  has  been  simplified,  recall  that  the  plow/soil  interaction  model  also  accounts  for 
the  fluid  phase  of  the  soil  and  remains  an  effective  stress  model. 

The  aforementioned  Drucker-Prager  vertex  singularity  problem  has  been  recognized  in 
the  literature  by  Abbo  and  Sloan  (1993).  Their  solution  is  to  approximate  the  classical  Drucker- 
Prager  yield  surface  with  a  smooth  hyperbola  that  is  continuous  and  differentiable  for  all  states  of 
stress  on  the  yield  surface  including  tensile  states.  The  hyperbola  can  be  fitted  to  the  Mohr- 
Coulomb  surface  (straight  part  of  the  Drucker-Prager  surface)  with  two  parameters.  The 
hyperbolic  approximation  everywhere  approaches  the  classic  Drucker-Prager  yield  surface  from 
the  inside,  and  so  the  approximation  of  soil  strength  is  conservative  with  this  method.  An 
approximate  Drucker-Prager  model  based  upon  this  comparatively  simple  fix  would  be  highly 
desirable  as  an  option  to  the  existing  Modified  Drucker-Prager  material  model,  since  no 
compromise  would  be  necessary  regarding  calibration  of  observed  strain  hardening  and  softening 
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behavior  in  sands.  (Such  a  fix  has  been  announced  by  HKS,  Inc.  for  ABAQUS/Standard, 
Version  5.5.) 

3.6.2  Tearing  the  Soil  in  Advance  of  the  Plow.  The  failure  mode  of  the  soil  directly  in 
front  of  the  plow  tip  is  an  important  modeling  consideration.  Two  failure  modes  were 
investigated  in  this  study;  a  displacement  mode,  and  a  tearing  mode  (Figure  3-4).  The 
displacement  failure  mode  is  comparatively  easy  to  simulate  and  was  implemented  first.  This 
mode  assumes  that  the  resistance  to  the  plow  through  the  soil  is  due  to  displacement  of  the  soil  to 
the  side,  and  that  in  the  direction  perpendicular  to  the  centerline,  the  x-direction,  only 
compression  stress  can  be  sustained  at  the  centerline.  It  is  implemented  via  a  unilateral  contact 
condition  in  the  x-direction  at  the  centerline,  allowing  centerline  nodes  to  freely  displace  in  the 
positive  x-direction  while  preventing  displacement  in  the  negative  x-direction,  in  advance  of  the 
plow  tip  (Figure  3-3).  This  is  implemented  with  an  ^INTERFACE  command  operating  on  a  set 
of  (specially  created)  centerline  elements  establishing  them  as  contact  elements,  and  with  a 
*  RIGID  SURFACE  option  establishing  a  rigid  surface  at  the  centerline  with  which  the  contact 
elements  are  to  interact. 

It  was  found  that  the  resulting  wedging  action  of  the  plow  tip  reduces  the  critical  surface 
over  which  the  soil  is  in  contact  with  the  blade  bevel  to  about  25  percent  of  the  total  surface.  In 
the  absence  of  friction,  plow  resistance  is  due  exclusively  to  bearing  stress  on  this  surface. 
Consequently,  prediction  of  this  resistance  was  very  low,  and  the  addition  of  side  friction 
between  the  soil  and  blade  would  be  the  only  way  to  correct  this  deficiency  associated  with  the 
displacement  failure  mode.  Further,  a  corresponding  prediction  of  a  very  distinct  gap  in  the  soil 
well  ahead  of  the  plow  tip  contradicts  experimental  observation.  Consequently,  this  failure  mode 
was  rejected,  even  though  it  was  not  apparent  at  the  time  how  the  alternative  tearing  mode  would 
be  modeled. 

The  tearing  failure  mode  assumes  that  the  plow  tip  is  the  main  source  of  resistance,  or 
conversely,  that  side  friction  is  a  minor  factor,  and  that  to  advance,  the  plow  must  fail  the  soil  by 
tearing  it  away  from  the  centerline.  This  tearing  failure  mode  is  also  known  in  the  study  of 
fracture  mechanics  of  brittle  materials  as  a  mode-I  failure.  Fortunately,  ABAQUS  provides  crack 
propagation  simulation  features  which  support  fracture  mechanics  studies,  and  moreover,  these 
features  are  also  available  in  conjunction  with  the  present  plow/soil  interaction  model.  However, 
this  additional  complexity  also  proved  to  degrade  the  performance  or  run-time  of  the  model. 

Adding  the  *BOND  SURFACE  sub-option  to  the  *  INTERFACE  modeling  option,  the 
prescribed  unilateral  contact  surface  along  the  centerline  will  then  also  provide  for  initially 
bonded  nodes.  These  nodes  may  also  debond,  as  prescribed  by  a  further  specification  of  a 
subsequent  *  DEBOND  load  history  option.  This  latter  option  defines  a  linear  rate  at  which  the 
stress  carried  by  the  bonded  interface  is  diminished  as  a  function  of  time  after  the  interface 
begins  to  debond.  Its  sub-option,  *CRACK  GROWTH,  is  used  to  specify  a  linear  crack  growth 
rate  beginning  at  a  fixed  reference  node,  CRACKREF .  Hence,  a  constant  crack  tip  velocity ,  Vc, 
with  which  the  interface  debonds  beginning  at  the  reference  node  is  specified.  The  prescribed 
velocity  should  be  only  slightly  greater  than  the  prescribed  velocity  of  the  plow,  V,  so  that  the 
plow  does  not  overtake  the  crack  tip,  and  the  interface  does  not  debond  too  far  in  advance  of  the 
plow,  forming  a  large  gap  as  described  for  the  displacement  failure  mode. 
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The  x-direction  force  at  a  bonded  node  prior  to  separation,  Fo,  is  reduced  linearly  to  zero 
over  an  interval  of  time,  ±to,  centered  at  the  time  the  crack  arrives  at  the  node.  The  time  tolerance 
to  is  user-prescribed.  In  this  study,  it  is  set  to  the  time  that  it  takes  the  plow  to  travel  half  an 
element  length,  Ay,  along  the  centerline.  Thus,  the  nodal  point  bond  force,  Fb,  reduces  to  zero  in 
an  interval  2to  centered  about  the  time  the  crack  arrives,  as  shown  in  Figure  3-11.  Values  for 
these  parameters  are  presented  in  Table  3-2,  with  the  crack  tip  velocity  which  is  selected  to  be 
less  than  0.1  percent  higher  than  the  prescribed  plow  velocity,  V. 

Table  3-2.  Tearing  Failure  Mode  Modeling  Parameters 


V 

Ay 

(in.) 

to 

(ms) 

Vc 

(in./sec) 

(fps) 

(in./sec) 

1 

12 

0.125 

12.01 

3 

36 

0.125 

-  ■  r 

ittBslii 

5 

60 

0.125 

1.040 

6 

72 

0.125 

0.868 

72.06 

8 

96 

0.125 

0.651 

96.08 

10 

120 

0.125 

0.521 

120.10 

13 

156 

0.125 

0.401 

156.13 

Though  Vc  is  only  slightly  greater  than  V,  the  crack  will  gradually  pull  away  from  the 

plow  tip  according  to  this  tearing  failure  model.  This  in  turn  causes  the  loaded  bearing  area  of  the 
plow  tip  to  diminish  slightly  with  time,  so  that  the  steady-state  resistance  will  correspondingly 
diminish  slightly  with  time.  Nonetheless,  the  tear  failure  model  accounts  for  substantially  larger 
(~2.5  times  larger)  predictions  for  plow  resistance  force  than  the  displacement  failure  model, 
primarily  because  it  coerces  a  larger  soil-bearing  area  at  the  plow  tip.  Since  the  Drucker-Prager 
model  is  pressure-dependent,  and  the  tearing  failure  model  also  tends  to  produce  higher  mean 
pressure  stress  at  the  plow  tip  (and  greater  negative  pore  pressure  as  well)  the  material  through 
which  the  plow  moves  is  correspondingly  stronger.  This  also  accounts  for  higher  predictions  of 
resistance. 

Despite  this  very  positive  additional  feature  of  the  plow/soil  interaction  model,  additional 
research  is  required  to  establish  the  actual  failure  characteristics  of  the  soil  at  the  plow  tip,  and 
then  to  model  the  actual  failure  mode. 

3.6.3  Modeling  Drainage  Paths  in  Two  Dimensions.  Plow  resistance  which  was 
dependent  on  plow  velocity,  as  observed  experimentally,  was  not  achieved  with  the  plow/soil 
interaction  model  as  described  above.  Further,  the  continuity  equation,  which  governs  the  pore 
fluid  phase,  and  which  the  model  obeys,  would  indicate  that  velocity  dependence  should  be  a 
characteristic  of  the  solution.  For  velocities  below  the  critical  velocity,  resistance  should  be  a 
function  of  the  rate  at  which  pore  fluid  flows  in  the  neighborhood  of  the  advancing  plow  tip. 

The  influence  of  the  continuity  equation  is  felt  only  when  the  fluid  is  given  a  chance  to 
flow.  If  the  plow  is  moving  so  fast  through  the  soil  that  the  pore  fluid,  as  modeled,  had 
insufficient  time  to  flow,  then  the  soil,  in  effect,  behaves  as  an  incompressible  medium.  This 
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appears  to  be  the  condition  in  the  soil  when  plow  velocity  exceeds  the  critical  velocity.  Under 
this  condition,  the  theoretical  result  is  that  resistance  would  be  independent  of  plow  velocity. 
This  type  of  behavior  characterized  the  results  with  the  initial  plow/soil  interaction  model 
described  above.  A  further  indication  of  this  was  that  initial  predictions  of  maximum  (negative) 
pore  pressures  were  two  orders  of  magnitude  larger  than  experimentally  observed  values  (recall 
that  ABAQUS  has  no  cavitation  cutoff  mechanism).  Moreover,  the  location  of  the  maximum 
negative  pressure  was  invariably  aft  of  the  plow  tip  and  acting  on  the  side  of  the  blade,  rather 
than  at  the  plow  tip  where  it  had  been  observed  experimentally  in  the  laboratory  apparatus  (see 
discussion  in  Kutter  and  Voss,  1994).  To  remedy  this,  attention  was  focused  on  compensating  for 
the  three-dimensional  drainage  paths  which  characterize  the  actual  or  physical  model  (Figure  1- 

There  are  two  differences  in  the  drainage  path  characteristics  between  the  numerical 
model  and  the  physical  model. 

1 .  The  plane  strain  numerical  model  drainage  paths  are  two-dimensional,  whereas  in  the 
physical  model  the  drainage  paths  are  three-dimensional  and  have  been  estimated  at  3 
to  6  inches  in  length. 

2.  The  numerical  model  is  confined  and  has  no  free  surface,  while  actual  flow  in  the 
neighborhood  of  the  plow  tip  occurs  in  the  presence  of  a  free  surface. 

To  account  for  these  differences,  two  adjustments  were  made  to  the  plow/soil  interaction 
model.  The  first  was  discussed  above  in  connection  with  Figure  3-3,  namely,  specification  of  a 
zero  pore  pressure  boundary  condition  along  the  boundary  line,  y  =  6  inches.  This  improved 
control  of  the  large  negative  pore  pressures,  but  it  is  not  sufficient  by  itself.  The  second 
adjustment  was  based  on  examination  of  some  experimental  soil  consolidation  results  by  Gibson 
et  al.  (1967).  In  this  work,  consolidation  tests  were  run  in  both  two-  and  three-dimensional  soil 
media,  and  consolidation  was  observed  to  occur  approximately  100  times  more  quickly  for  the 
three-dimensional  medium.  To  replicate  this  difference  in  the  plow/soil  interaction  model,  values 
for  soil  permeability  were  increased  two  orders  of  magnitude  beyond  experimentally  observed 
values  for  the  soil  in  the  physical  model. 

The  actual  coefficient  of  permeability  k  for  the  soil  was  empirically  determined  to  be  a 
function  of  void  ratio  cubed  (True,  1994), 

k  =  0.087  e^  (in/sec) 

as  shown  in  Figure  3-12  for  the  anticipated  range  of  void  ratio.  The  prescribed  value  of 
permeability  for  the  plow/soil  interaction  model  was  given  by  this  relation,  times  a  factor  of  100, 
to  account  for  drainage  path  differences  in  two-  and  three-dimensional  soil  consolidation  media 
according  to  the  aforementioned  experimental  result. 

Results  from  the  model  subsequent  to  this  change  were  immediately  recognizable  and 
positive.  The  maximum  predicted  (negative)  pore  pressures  dropped  by  two  orders  of  magnitude 
to  levels  representative  of  the  laboratory  test  apparatus,  and  the  maximum  pore  pressure 
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distribution  enveloped  the  plow  tip,  consistent  with  observed  pore  pressure  distribution  in  the 
physical  model. 

3.6.4  Solution  Control  Parameters.  A  summary  of  the  features  and  size  of  the 
numerical  plow/soil  interaction  model  is  given  in  Table  3-3.  In  addition  to  the  usage  of  a 
Drucker-Prager  material  model,  geometric  and  material  nonlinearity,  coupled 
deformation/porous  flow  behavior,  and  interface  nonlinear  behavior,  the  model  has  8,787 
unknowns  (nodal  displacements  and  pore  pressures).  With  this  complexity,  solution  convergence 
problems  may  well  be  anticipated.  It  is  then  useful  to  relax  certain  solution  control  parameters. 
ABAQUS  provides  for  this  with  a  *CONTROL  command.  Further,  default  tolerance  values  are 
considered  overly  tight,  and  often  engineering  accuracy  can  be  retained  after  these  tolerances  are 
relaxed  to  foster  improved  performance.  In  the  present  study  it  was  necessary  to  use  non-default 
values  of  residual  and  correction  quantities  for  both  the  solid  and  fluid  phases. 

The  default  and  non-default  solution  control  parameters  which  were  used  are  given  in 
Table  3-4.  Default  values  for  residual  force  and  displacement  correction,  5.0E-03  and  0.01, 
respectively,  have  been  overridden  to  promote  solution  convergence  in  the  solid  phase.  The 
largest  residual  nodal  point  force  in  any  solution  increment  must  be  no  larger  than  0.01  or  1 
percent  of  the  largest  increment  in  nodal  point  force.  Further,  the  ratio  of  largest  solution 
correction  (for  nodal  displacement)  to  largest  solution  increment  (for  nodal  displacement)  must 
be  less  than  0.5.  Otherwise,  convergence  is  not  achieved  for  the  increment  and  further  iteration 
occurs.  All  other  solution  control  parameters  were  assigned  default  values. 
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Figure  3-11.  Centerline  nodal  point  debond  force. 
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Figure  3-12.  Permeability  variation  with  void  ratio. 


Table  3-3.  Summary  of  Finite  Element  Features  for  2-D  Plow/Soil  Interaction  Model 


Type  of  soil  element 

8-node  plane  strain  quadrilateral  element  w/bi- 
quadratic  displacement  and  bilinear  pore  pressure 
interpolation. 

No.  of  soil  elements 

1,274 

No.  of  nodal  points 

3,935 

No.  of  rigid  surfaces 

2;  Rigid  wall  centerline  and  rigid  blade. 

No.  of  propagating  cracks 

1 ;  Constant  velocity  crack  forming  at  centerline  ahead 
of  plow  tip  with  linearly  debonding  nodal  force. 

Type  of  contact  element 

2-node  planar  element  w/linear  interpolation. 

No.  of  contact  elements 

149;  74  on  wall  centerline  and  75  on  soil  centerline. 

Type  of  multipoint  constraint 

Lagrange  multiplier;  linear  interpolation  of  soil 
element  midside  node  displacement  along  centerline. 

No.  of  multipoint  constraints 

75 

Pore  pressure  boundary  constraints 

Pore  pressure  constrained  to  zero  at  151  centerline 
nodes,  and  at  151  edge  nodes  (line,  x  =  6  inches). 

Displacement  boundary  constraints 

Normal  displacement  constrained  to  zero  at  center-line, 
and  along  edge  lines  x  =  6  inches  and  y  =  12  inches. 

No.  of  degrees  of  freedom  (DOF) 
Maximum  DOF  wavefront  estimate 

8,787  (including  Lagrange  Multipliers) 

107  (optimized;  RMS  estimate,  104) _ 

Table  3-4.  Default  and  Non-Default  Convergence  Parameter  Settings  for  Plow/Soil 

Interaction  Numerical  Model 


Parameter 

Value 

Convergence  Tolerance  Parameters  for  Force 

Criterion  for  residual  force  for  a  nonlinear  problem3 

Criterion  for  displacement  correction  in  a  nonlinear  problem* 
Initial  value  of  time  average  force 

Alternate  crit.  for  residual  force  for  a  nonlinear  problem 

Criterion  for  zero  force  relative  to  time  average  force 

1.000e-02 

5.000e-01 

1.000e-02 

2.000e-02 

1.000e-05 

Convergence  Tolerance  Parameters  for  Volume  Flux 

Criterion  for  residual  vol.  flux  for  a  nonlinear  problem* 

Criterion  for  p.press.  correction  in  a  nonlinear  problem* 

Initial  value  of  time  average  volume  flux 

Alternate  crit.  for  residual  vol.  flux  for  a  nonlinear  prob. 

Criterion  for  zero  vol.  flux  relative  to  time  avrg.  vol.  flux 

5.000e-03 

1.000e-02 

1.000e-02 

2.000e-02 

1.000e-05 

a  Indicates  use  of  non-default  convergence  tolerance  parameter  selected  with  a  *  CONTROLS  command 
(e.g.,  *  CONTROLS,  PARAMETERS=FIELD,  FIELD=DISPLACEMENT  with  the  following  data  card; 
0.0 1”  0.5). 


52 


4.0  COMPUTATIONAL  RESULTS  OF  PLOW/SOIL  INTERACTION  MODEL 


To  gain  an  understanding  of  the  multistep  numerical  model  of  plow/soil  interaction, 
typical  quasi-static  results  are  presented  and  discussed  in  some  detail.  Results  are  presented 
primarily  in  graphical  form  using  color  to  promote  their  interpretation.  These  results  include  unit 
plow  resistance  forces  and  deformation,  stress  and  pore  pressure  in  the  surrounding  soil  during 
both  transient  and  steady-state  motion  conditions.  Time  histories  and  spatially  contoured 
distributions  are  presented  for  most  response  variables.  In  the  former  case,  data  are  reported  for  a 
series  of  node  points  and  elements  along  the  centerline  or  path  of  travel  for  the  plow.  In  the  latter 
case,  data  are  reported  for  both  the  solid  and  fluid  phases  of  soil  surrounding  the  tip  of  the  plow. 
These  data  are  reported  from  high  order  interpolated  elements  for  the  solid  phase,  and  from  low 
order  interpolated  elements  for  the  pore  fluid  phase.  They  typify  the  solution  detail  of  which  the 
finite  element  method  is  capable  for  complex  nonlinear  problems.  The  engineering  parameter 
studies  derived  from  the  model  and  reported  in  Section  5.0  are  based  on  these  data. 

4.1  Model  Displacements 

It  is  useful  to  first  study  the  soil  displacement  of  the  plow/soil  interaction  model  to 
become  familiar  with  its  behavior.  Also,  soil  deformation  is  directly  related  to  how  a  trench  is 
formed,  which  is  of  practical  interest.  In  Figure  4-1,  the  displacement  of  the  plow  blade  (6.5 
inches  in  length)  and  the  surrounding  soil  are  shown  when  the  plow  has  penetrated  into  the  soil. 
Corresponding  response  time  histories  of  resisting  force  verify  that  the  plow  is  moving  under 
steady-state  conditions  when  having  penetrated  this  distance.  The  characteristics  of  soil 
deformation  shown  are  typical  of  the  plow/soil  interaction  model,  irrespective  of  the  steady-state 
plow  velocity.  The  bevel  of  the  plow  tip  bears  directly  against  the  soil  pushing  it  aside  as  the 
plow  blade  progresses.  The  deformation  in  the  soil  medium  is  maximum  in  the  region  of  the 
plow  tip.  Soon  after  the  shoulder  of  the  plow  tip  moves  past,  the  soil  separates  from  the  plow 
blade.  Soil  along  the  plow  blade  tends  to  shear  in  the  direction  of  travel.  Due  to  local  plastic 
deformation  of  the  soil,  interface  separation  occurs  and  appears  to  increase  with  distance  along 
the  blade  to  a  maximum  at  the  aft  end.  A  distinct  wake  is  formed  aft  of  the  plow  which  is  also  in 
qualitative  agreement  with  experimental  observation.  Rebound  is  evident  in  the  wake  indicating 
elastic  deformation  in  the  soil  mass  as  well.  Aft  of  the  plow,  a  substantial  amount  of  extruded 
soil  material  is  shown  which  is  actually  flowing  opposite  to  the  plow  in  the  present  plow/soil 
interaction  model.  This  material  is  escaping  the  confinement  of  the  model’s  rigid  boundaries  and 
is  displaced  by  the  penetrating  plow  blade.  It  is  an  artifact  of  the  model. 

A  section  of  the  model  along  the  path  of  travel  is  presented  in  Figure  4-2  to  show  the 
location  of  nodes  and  elements  where  response  time  histories  are  sampled  to  support  discussion 
of  response  of  the  system.  Data  are  sampled  at  the  node  points  and  the  elements  shown  which  are 
located  directly  on  or  adjacent  to  the  centerline  of  the  path  of  the  plow.  It  will  be  useful  to  refer 
to  these  locations  when  studying  strain,  stress  and  pore  pressure  results.  Up  to  1,000  implicit 
solution  time  increments  are  often  required  to  complete  a  typical  analysis.  Data  is  written  to 
output  files  every  10  increments,  so  time-history  graphs  are  often  constructed  with  up  to  100  data 
points. 

In  the  following  discussion,  the  plow  is  moving  at  a  steady-state  velocity  of  5  fps,  the 
blade  (full)  width  is  0.5  inches,  the  blade  length  is  6.5  inches,  the  included  angle  of  the  beveled 
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tip  is  60  degrees  (Figure  1-lb),  and  the  plow/soil  interface  is  frictionless.  Depth  and  sweep  angle 
of  the  blade  leading  edge  are  not  parameters  of  this  plane  strain  plow/soil  interaction  model,  and 
are  therefore  absent  from  the  discussion.  The  depth  is  assumed  to  be  unity  and  the  sweep  angle  is 
assumed  to  be  90  degrees. 

In  summary,  the  computed  soil  displacements  evidence  the  basic  behavior  of  the 
plow/soil  interaction  model.  The  important  displacements  are  confined  to  the  plow  tip  region. 
Essentially,  only  the  plow  tip  is  in  contact  with  the  soil,  and  the  response  of  the  model  is 
otherwise  independent  of  the  length  of  the  plow  blade.  Specific  features  of  model  displacements 
are: 


•  Maximum  localized  deformation  occurs  in  the  region  of  the  plow  tip. 

•  Separation  occurs  at  a  point  just  aft  the  shoulder  of  the  plow  tip,  and  increases  with 
distance  along  the  plow  blade. 

•  Soil  is  sheared  in  the  direction  of  motion  and  contains  both  elastic  and  plastic 
deformation. 

•  A  permanent  wake  is  formed  behind  the  plow. 

4.2  Plow  Forces 

Since  the  plow/soil  interaction  model  is  not  three-dimensional,  the  total  steady-state 
towing  resistance  cannot  be  computed  directly  and  would  have  to  be  estimated.  The  estimation 
method  is  crude  and  includes  the  computed  force  history  per  unit  depth  given  by  model,  a 
reasonable  interpolation  of  bearing  stress  through  the  depth  of  the  soil,  and  experimental 
observation  of  the  tow  force.  These  considerations  are  discussed  along  with  the  computed  contact 
stress  distribution  on  the  plow  blade. 

4.2.1  Unit  Tow  Force  History.  Time  histories  of  the  two  rectangular  components  of 
force-per-unit  depth  acting  on  the  plow  are  shown  in  Figure  4-3.  The  computed  time  history  of 
the  force-per-unit  depth  required  to  tow  the  (half)  blade  through  the  soil  at  5  fps  is  shown  in 
Figure  4-3a.  (Its  magnitude  can  be  corrected  for  the  full  blade  thickness  by  multiplying  by  two.) 
Because  the  two-dimensional  plow/soil  interaction  model  has  unit  thickness  in  the  depth 
direction,  the  tow  force  is  a  force  per-unit-depth  having  dimensions,  lb/in.  It  is  the  integral  of  in¬ 
line  bearing  stress  (total  stress)  over  the  projected  edge  of  the  blade.  In  these  data  the  initial  or 
transient  phase  of  the  tow  force  is  clearly  evident  and  shows  a  slight  overshoot  involving  initial 
penetration  of  the  plow  blade,  after  which  a  steady-state  force  condition  occurs.  The  oscillations 
in  the  data  are  generally  due  to  numerical  noise  relating  to  the  numerical  algorithms  for  handling 
interface  contact  and  nonlinear  material  behavior.  (They  are  not  unlike  the  noise  oscillations  seen 
in  typical  graphs  of  measured  force  time  histories  in  the  laboratory.)  The  steady-state  magnitude 
of  the  unit  tow  force  can  be  obtained  by  merely  averaging  these  data  near  the  end  of  the  range 
shown  when  the  entire  blade  length  is  sufficiently  embedded  (Figure  4-1). 
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Figure  4-1 .  Displacement  of  plow  and  surrounding  soil  elements. 
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(b)  Lateral  and  tow  (half)  force  components  per  unit  depth  of  blade. 


4.2.2  Unit  Lateral  Plow  Force  History.  In  addition  to  the  tow  force  or  in-line 
component,  the  plow/soil  interaction  model  computes  the  unit  lateral  or  side  component  of  force 
on  the  plow  which  acts  perpendicularly  to  the  blade  (tending  to  pinch  the  blade  )  as  shown  in 
Figure  4-3b.  This  component  is  the  integral  of  total  soil  stress  over  the  area  of  the  leading  edge  of 
the  blade  up  to  the  point  of  interface  separation  (Figure  4-1).  Graphed  against  the  half-tow  force 
component,  the  lateral  force  is  approximately  twice  as  large,  which  makes  the  lateral  and  the  tow 
force  components  about  equal.  Both  components  clearly  evidence  steady-state  conditions. 
Depending  on  blade  roughness  and  interface  conditions,  the  lateral  component  could  lead  to 
significant  side  friction  force.  For  the  frictionless  blade  condition,  this  component  has  no  role  in 
the  estimate  of  resistance  force.  Instrumentation  in  the  laboratory  apparatus  cannot  measure  the 
lateral  component;  it  can  only  measure  the  unbalanced  force  in  the  lateral  direction. 

4.2.3  Unit  and  Total  Plow  Resistance.  A  total  resistance  force  for  the  plow,  F,  can  be 
crudely  estimated  in  conjunction  with  the  steady-state  unit  tow  force  component  computed  by  the 
plow/soil  interaction  model  (Figure  4-3),  herein  denoted  fSs-  However,  the  main  points  of  this 
exercise  are  to  emphasize  that  the  two-dimensional  plow/soil  interaction  model  cannot  predict 
total  resistance  force,  and  to  reinforce  the  meaning  of  the  quantity,  unit  plow  force.  Two 
assumptions  are  necessary,  as  depicted  in  Figure  4-4: 

1 .  The  plow/soil  interaction  model  is  assumed  to  pertain  to  a  two-dimensional  region  of 

soil  at  an  arbitrary  depth,  H’,  between  the  surface  of  the  soil  bed  and  the  bottom  of 
the  plow  at  depth  H.  The  value  of  H’  may  be  determined  or  calibrated  in  conjunction 
with  measured  data  for  total  resistance  force,  Fr?  so  as  to  achieve  agreement 
between  measured  and  estimated  resistance  force. 

2.  The  distribution  of  bearing  or  resistance  force  with  depth  along  the  leading  edge  of  the 
plow  cannot  be  predicted  by  a  two-dimensional  plow/soil  interaction  model,  and  must 
be  assumed.  Experimental  evidence  from  the  laboratory  apparatus  suggests  that  the 
total  resistance  force  acts  below  the  center  of  the  plow,  and  erosion/wear  of  the 
leading  edge  of  blades  tested  in  the  field  shows  an  increase  with  depth  along  the 
leading  edge.  Thus,  for  convenience,  a  linear  distribution  of  bearing  force  with  depth 
may  be  assumed. 

In  the  present  case  of  a  frictionless  1 /2-inch  blade  with  depth  H  =  6  inches  (Figure  1-lb) 
and  moving  at  5  fps  through  a  100  percent  saturated  sand,  the  computed  unit  steady-state  tow 
force  fss  is  approximately  126  lb/in.  (Figure  4-3  a),  and  the  measured  total  resistance  force  from 
the  laboratory  apparatus,  Fr?  is  approximately  950  pounds  (Figure  l-2a).  The  corresponding 
calculated  value  of  H’  would  therefore  be  about  2.4  inches,  using  the  formulas  in  Figure  4-4. 
This  value  can  then  be  considered  a  constant,  and  knowing  fss  from  the  model,  a  total  resistance 
force  F  can  then  be  estimated,  albeit  very  crudely. 
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4.2.4  Lateral  Contact  Force  Distribution.  Lateral  contact  force  is  a  force-per-unit 
length  along  the  plow/soil  interface,  and  it  is  based  on  total  soil  stress.  The  lateral  contact  force 
distribution  over  the  plow  blade  is  graphed  in  Figure  4-5.  These  data  are  reported  by  the  contact 
elements  located  at  the  plow/soil  interface.  The  lateral  force-per-unit  depth  acting  on  the  plow, 
which  is  independently  reported  by  the  rigid  plow  blade  element  (Figure  4-3b),  is  the  integral 
over  the  contact  zone  of  the  contact  force  distribution  shown.  Profiles  of  both  the  deformed  soil 
and  rigid  blade  surfaces  are  also  shown.  Axes  with  tick  marks  have  their  origins  along  the 
deformed  soil  interface.  This  interface  is  coincident  with  the  rigid  plow  blade  for  only  a  short 
distance  near  the  plow  tip  where  contact  exists.  The  contact  force  is  zero  along  the  interface  aft  of 
the  blade  shoulder  (and  forward  of  the  blade  tip)  and  non-zero  near  the  tip  of  the  plow.  The  peak 
lateral  contact  force  occurs  at  the  shoulder  of  the  blade  tip.  Magnitudes  shown  are  indicative  of 
steady-state  conditions,  but  they  are  greater  under  transient  or  startup  conditions,  just  as  plow 
forces  are  greater  under  transient  conditions  (Figure  4-3). 

4.2.5  Summary  of  Model  Plow  Forces.  Though  it  can  be  measured  adequately  in  the 
laboratory  apparatus,  it  is  not  possible  with  a  two-dimensional  plane  strain  model  of  plow/soil 
interaction  to  directly  predict  the  total  resistance  force  acting  on  the  plow.  However,  plow  forces 
per  unit  depth,  or  unit  plow  forces,  can  be  predicted  which  contribute  to  the  understanding  of 
plow/soil  interaction,  as  follows: 

•  It  is  more  accurate  to  speak  in  terms  of  plow  tip  forces,  rather  than  plow  forces,  for  all 
mechanical  interaction  with  the  soil  corresponds  to  the  plow  tip. 

•  With  the  plow  undergoing  constant  velocity  motion,  time  histories  of  plow  tip  forces 
exhibit  two  phases  of  response,  initial  startup  transient  forces  as  the  blade  enters  the 
soil,  and  subsequent  steady-state  forces  after  the  blade  is  fully  embedded.  The  peak 
transient  force  overshoots  the  steady-state  force,  and  the  model  predicts  their  relative 
magnitudes. 

•  The  two  rectangular  components  of  plow  tip  force  are  predicted;  i.e.,  the  in-line  or 
tow  force  component  and  the  lateral  or  side  force  component.  The  model  predicts 
their  relative  difference,  whereas  the  lateral  force  component  cannot  be  measured  in 
the  current  laboratory  apparatus. 

•  Lateral  contact  forces  along  the  bevel  of  the  plow  tip  vary  approximately  linearly 
from  the  tip  to  a  maximum  value  at  the  shoulder. 
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Figure  4-4.  Method  for  extrapolating  to  total  resistance  force  on  plow. 
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Figure  4.5.  Distribution  of  lateral  interface  contact  force  on  plow  at  5  fps 
steady-state  plow  velocity. 


4.3  Deformation  in  the  Soil 


The  plasticity  constitutive  model  (Section  3.5)  employed  for  modeling  the  soil  is  an 
important  aspect  of  the  plow/soil  interaction  model.  Investigating  the  computed  deformation  in 
the  soil  is  also  important  because,  among  other  things,  the  soil  is  a  dilatant  material,  and  the  soil 
model  should  evidence  this  characteristic.  To  analyze  the  computed  volume  increase  that  occurs 
as  the  soil  is  sheared,  response  of  the  direct  or  normal  strain  component  must  be  discussed.  For 
completeness,  the  shear  strains  are  also  included  in  the  discussion.  In  classic  plasticity  theory, 
strain  is  resolved  into  an  elastic  and  a  plastic  component.  The  following  data  will  show  that,  near 
the  plow/soil  interface,  the  plastic  strain  component  dominates  the  elastic  strain  component  for 
both  the  direct  and  shear  components. 

Strains  are  element  data  and  the  data  presented  are  computed  at  the  centroid  of  each 
element.  In  the  former,  the  data  are  taken  from  the  selected  centerline  elements  (Figure  4-2),  and 
in  the  latter,  the  strain  fields  shown  are  enlarged  views  of  a  region  near  the  plow  tip  when  the 
plow  is  at  a  position  undergoing  steady-state  motion  (Figure  4-1).  Regarding  the  strain  field 
contour  data,  it  is  important  to  point  out  that  the  post-processor  algorithm  averages  the  element 
centroid  values  at  adjacent  node  points  before  constructing  the  contour  plots. 

4.3.1  Direct  Strain.  The  two  non-zero  components  of  direct  strain,  Ell  and  E22,  are 
presented  in  Figure  4-6,  and  each  of  these  strain  components  is  the  sum  of  an  elastic  and  a  plastic 
part.  The  third  component,  E33,  is  constrained  to  zero  under  the  assumption  of  a  plane  strain 
model.  The  total  strain  and  plastic  strain  response  are  presented;  elastic  strain  can  be  inferred  as 
their  difference. 

Time  histories  for  Ell,  the  lateral  component  of  total  strain,  along  the  centerline  (CL)  or 
path  of  the  plow,  are  shown  in  Figure  4-6a.  This  strain  component  is  apparently  tensile  along  the 
centerline.  The  history  traces  are  characterized  by  a  peak  response  which  occurs  during  arrival 
and  departure  of  the  plow  tip,  and  thereafter  subside  to  a  steady  residual  tensile  strain  response 
which  is  a  substantial  fraction  of  the  peak  response.  The  residual  strain  characterizes  the  total 
lateral  strain  component  in  the  soil  at  the  free  edge  of  the  groove  or  wake  produced  by  the  plow 
tip.  The  spatial  distribution  of  this  strain  under  steady-state  conditions,  shown  in  Figure  4-6b, 
indicates  that  actual  peak  strains  occur  just  off  the  interface  and  adjacent  to  the  plow  tip  rather 
than  on  the  interface.  Residual  strain  extends  to  a  distance  of  about  two  blade  thickness  to  the 
side  of  the  groove,  beyond  which  it  quickly  dissipates. 

Time  histories  for  PE11,  the  plastic  part  of  the  lateral  strain  component  along  the 
centerline,  are  shown  in  Figure  4-6c.  Compared  to  the  total  strain  (Figure  4-6a),  the  plastic  strain 
is  nearly  equal  in  the  residual  portion  of  the  response.  Thus,  the  residual  elastic  strain  is  minute 
in  these  data,  as  would  be  expected  since  the  soil  has  unloaded.  The  spatial  distribution  of  PEI  1 
which  is  shown  in  Figure  4-6d  is  also  very  much  like  that  for  El  1  (Figure  4-6b).  Most  of  the  total 
strain  component  Ell  is  plastic  strain  and,  therefore,  there  should  be  little  rebound  in  the  soil  at 
the  edge  or  lip  of  the  soil  groove  after  the  blade  shoulder  passes. 

Time  histories  for  E22,  the  component  of  total  strain  in  the  direction  of  the  plow  motion 
along  the  centerline  or  path  of  the  plow,  are  shown  in  Figure  4-6e.  The  history  traces  are 
characterized  by  a  transient  compressive  strain  response  which  occurs  during  arrival  and 
departure  of  the  plow  tip,  and  thereafter  reverses  to  a  steady  residual  tensile  strain  response 
which  is  comparable  in  magnitude  to  the  transient  response. 
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Figure  4-6.  Direct  strain  components  in  soil  at  5  fps  steady-state  plow  velocity. 
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(b)  Total  strain  El  1  distribution  near  plow  tip. 


LIME  VARIABLE 
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(c)  Plastic  strain  PEI  1  histories  at  sampled  response  points  along  path  of  plow. 
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(d)  Plastic  strain  PEI  1  distribution  near  plow  tip. 


Both  direct  strain  components  at  the  edge  of  the  groove  are  tensile  in  the  residual 
response  and,  therefore,  the  material  has  dilated  due  to  the  shearing  action  of  the  plow.  This 
would  be  consistent  with  a  dilatant  soil.  The  spatial  distribution  of  E22  under  steady-state 
conditions,  shown  in  Figure  4-6f,  indicates  that  actual  maximum  values  occur  just  off  the 
interface  adjacent  to  the  plow  tip  rather  than  on  the  interface,  similar  to  El  1.  Residual  strain  E22 
is  slightly  larger  than  Ell  along  the  groove.  However,  maximum  transient  values  for  Ell  are 
larger  than  for  E22.  Except  for  the  plastic  region  near  the  plow  tip  and  also  along  the  groove, 
both  direct  strain  components  are  negligible  in  the  surrounding  soil.. 

Time  histories  for  PE22,  the  in-line  component  of  plastic  strain  along  the  centerline,  are 
shown  in  Figure  4-6g.  When  compared  to  the  total  strain  E22  (Figure  4-6e),  the  plastic 
component  is  also  compressive,  behaves  in  nearly  the  same  manner,  and  also  accounts  for  most 
of  the  total  strain.  The  spatial  distribution  of  PE22  which  is  shown  in  Figure  4-6h  is  also  very 
much  like  that  for  E22  (Figure  4-6f).  Most  of  both  residual  direct  strain  components  El  1  and  E22 
are  plastic  along  the  interface. 

4.3.2  Volume  Change.  Volume  change  in  the  material  is  reported  in  terms  of  time 
histories  of  element  volume  for  the  selected  elements,  as  shown  in  Figure  4-7.  Each  of  these 
elements  are  0.125  inches  square  in  the  undeformed  configuration,  with  a  constant  unit  depth. 
Therefore,  their  initial  or  undeformed  volume  is  0.0156  inches  cubed.  Element  volume  increases 
with  arrival  of  the  plow  blade.  During  arrival  and  departure  of  the  plow  tip,  there  is  a  distinct 
initial  peak  increase  with  a  small  rebound,  and  then,  as  the  blade  shoulder  slides  by  shearing  the 
soil,  the  element  volume  rises  to  a  maximum  residual  level.  In  these  data,  the  residual  volume 
expansion  is  about  3 1  percent  and  is  due  to  plastic  tensile  strain. 

Alternatively,  if  element  volume  is  denoted  as  EVOL,  then  the  change  in  element 
volume,  AEVOL,  equals  the  algebraic  sum  of  the  direct  strains  El  1  and  E22.  For  the  residual 
response,  this  sum  is  about  27  percent  (Figures  4-6a  and  4-6e).  The  average  strain  is  therefore 
reasonably  within  the  calibration  range  for  the  dilation  angle  (Figure  3.9b)  of  the  Drucker-Prager 
model. 


4.3.3  Shear  Strain.  Time  histories  of  the  shear  strain,  E12,  along  the  centerline,  shown 
in  Figure  4-8(a),  exhibit  a  reversal  response  during  arrival  and  departure  of  the  plow  tip,  first 
positive  and  then  negative.  Due  to  symmetry,  the  shear  strain  must  remain  zero  at  the  centerline 
until  the  plow  tip  arrives  and  cleaves  the  soil,  and  then  the  response  includes  a  precipitous 
reversal  of  shear  strain.  The  residual  response  subsequent  to  passing  of  the  blade  shoulder  is 
negative  and  nearly  equal  to  the  peak  negative  response.  The  spatial  distribution  of  the  steady 
response  of  El 2  is  shown  in  Figure  4-8b.  Shear  strain  is  evident  in  the  distortion  of  the  element 
grid  lines  shown  which  otherwise  intersect  at  right  angles  in  the  undeformed  state.  The  actual 
peak  or  maximum  response  occurs  well  off  the  interface,  perhaps  two  blade  widths  away,  in  the 
region  of  the  plow  tip.  It  is  also  opposite  in  sign  to  the  shear  strain  which  is  located  in  an 
apparent  thin  interface  boundary  layer  of  soil  adjacent  to  the  plow  blade. 
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(e)  Total  strain  E22  histories  at  sampled  response  points  along  path  of  plow. 
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(h)  Plastic  strain  PE22  distribution  near  plow  tip. 
Figure  4-6.  Continued. 
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(b)  Total  shear  strain  El 2  distribution  near  plow  tip. 
Figure  4-8.  Continued. 


The  plastic  component  of  shear  strain,  PE  12,  is  nearly  identical  to  the  total  shear  strain  as 
shown  in  Figures  4-8c  and  d.  Thus,  the  elastic  shear  strain  component  is  essentially  zero  in  these 
data.  Plastic  strain  completely  characterizes  the  response  of  the  soil,  within  one  or  two  blade 
widths  from  the  blade.  Plastic  shear  strain,  PE  12,  dominates  residual  strain  in  the  wake  or  groove 
material,  while  direct  plastic  strain  dominates  in  the  region  nearest  the  plow  tip. 

4.3.4  Summary  of  Model  Strains  in  the  Soil.  Results  for  model  strain  reveal  a  residual, 
plastic  strain  response  for  the  soil  bordering  the  wake  formed  behind  the  plow.  Time  histories 
show  that  soil  develops  peak  strain  upon  arrival  of  the  blade  shoulder,  and  retains  a  sizable 
percentage  as  residual  strain  after  it  passes.  The  rise  and  fall  in  strain  is  precipitous.  Primary  soil 
strain  response  is  confined  to  one  or  two  blade  widths  from  the  plow,  and  strain  gradients  are 
large.  Specific  results  are: 

•  The  lateral  direct  strain  component  remains  tensile,  while  the  in-line  direct  strain 
component  is  compressive  forward  of  the  plow  tip,  and  reverses  to  a  tensile  residual 
response  aft  the  plow  tip.  Shear  strain  also  reverses  as  the  plow  tip  passes,  and  also 
retains  a  residual  response. 

•  Both  direct  strain  components  retain  a  significant  tensile  residual  response  in  the 
wake,  which  demonstrates  the  model  correctly  replicates  volumetric  expansion 
behavior  in  a  dilatant  soil  subjected  to  the  shearing  action  of  the  plow. 

4.4  Effective  Stress  in  the  Soil 

Direct  effective  stress  in  the  soil  is  the  algebraic  sum  of  the  total  direct  stress  and  the  pore 
pressure,  while  effective  and  total  shear  stresses  are  identical  quantities.  Under  plane  strain 
conditions  there  are  four  non-zero  effective  stress  components,  SI  1,  S22,  S33,  and  S12  which 
define  the  state  of  stress  in  the  soil  skeleton  at  any  point.  The  first  three  are  direct  stress 
components  and  the  fourth  is  the  shear  stress  component.  SI  1  and  S22  are  bearing  stresses,  but 
only  the  S22  component  contributes  directly  to  plow  resistance  force  since  friction  is  neglected  in 
the  present  analysis.  The  S33  component  arises  primarily  because  of  the  constraint  on 
deformation  inherent  in  the  otherwise  artificial  plane  strain  conditions. 

Two  derived  stresses  or  equivalent  stresses  are  also  reported.  Both  are  invariant 
quantities.  The  Mises  equivalent  stress  is  a  measure  of  deviatoric  (distortion)  behavior  in  the  soil 
skeleton.  The  equivalent  pressure  stress,  is  indicative  of  the  dilatation  behavior  in  the  soil 
skeleton.  In  this  study,  it  is  defined  as  the  negative  of  mean  effective  stress  or  hydrostatic  stress. 

Like  strains,  stresses  are  element  data,  and  the  data  presented  are  computed  at  the 
centroid  of  an  element,  whether  they  are  graphs  of  time-histories  of  stress  or  contour  plots  of 
spatial  distributions  of  stress.  In  the  former,  the  data  are  reported  at  selected  centerline  elements 
(Figure  4-2)  and  in  the  latter,  the  stress  fields  shown  are  enlarged  views  of  a  region  near  the  plow 
tip  when  the  plow  is  in  steady-state  motion  (at  nearly  the  position  shown  in  Figure  4-1).  In  the 
stress  contour  data,  the  post-processor  algorithm  averages  the  computed  element  centroid  values 
at  each  adjacent  node  point  before  constructing  the  contours. 
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The  important  peak  or  maximum  response  in  stress  history  data  is  associated  with  arrival 
and  departure  of  the  blade  shoulder  as  for  the  corresponding  strain  data.  Residual  stress  in  the 
wake  of  the  plow,  unlike  residual  strain,  is  invariably  zero. 

4.4.1  Direct  Stress.  The  results  for  direct  stress  components  SI  1  and  S22  are  shown  in 
Figure  4-9.  The  lateral  component  SI  1  is  compressive  along  the  centerline  or  path  of  travel  of  the 
plow  blade  according  to  the  time  histories  shown  in  Figure  4-9a.  These  graphs  rise  and  fall 
exactly  according  to  the  arrival  and  departure  of  the  tip  of  the  plow  blade.  Steady-state  behavior 
is  evident  from  these  data  at  all  centerline  elements  beyond  element  40.  Their  spatial  distribution 
in  Figure  4-9b  also  includes  a  graph  of  Sll  along  the  centerline  which  is  essentially  a  normal 
section  through  the  contour  map  along  the  interface.  This  graph  includes  a  line  drawn  through  the 
grid  axes  just  above  their  midpoint  that  indicates  the  zero-stress  datum.  The  graph  shows  that 
S 1 1  is  maximum  at  the  shoulder  of  the  blade  and  tapers  linearly  towards  zero  at  the  tip  of  the 
blade,  and  that  it  is  zero  forward  of  the  tip  along  the  centerline  and  also  aft  of  the  point  of 
interface  separation.  These  data  are  directly  related  to  the  aforementioned  contact  force  (Figure 
4-5)  where  the  approximately  linear  variation  over  the  bevel  of  the  blade  tip  is  perhaps  seen  more 
clearly.  This  lateral  stress  component  could  be  important  to  both  the  resistance  force  of  the  plow 
and  wear  of  the  blade  shoulder,  depending  on  the  coefficient  of  friction  (assumed  zero  here) 
between  the  plow  and  the  soil  (See  Section  5.2). 

The  in-line  component  S22  contributes  directly  to  the  unit  resistance  force.  Time  history 
graphs  are  shown  in  Figure  4-9c.  They  show  that  steady-state  behavior  is  achieved  at  all 
centerline  elements  beyond  element  50.  These  graphs  show  a  gradual  increase  in  S22  before  the 
blade  tip  arrives,  unlike  the  more  precipitous  graph  for  S 1 1 .  When  the  blade  tip  does  arrive  and 
depart  the  corresponding  increase  and  decrease  in  S22  are  also  precipitous.  The  peak  magnitude 
of  S22  is  noticeably  smaller  than  for  Sll.  These  data  (Sll  and  S22)  also  show  that  the  peak 
startup  or  transient  stresses  are  substantially  higher  than  the  maximum  steady-state  stresses.  In 
Figure  4-9(d)  the  spatial  distribution  of  S22  shows  only  a  minor  concentration  at  the  shoulder  of 
the  blade,  i.e.,  a  more  uniform  distribution  along  the  bevel.  Overall,  the  maximum  value  for  S22 
is  located  on  the  bevel  surface,  and  is  somewhat  smaller  than  the  maximum  value  for  Sll.  The 
predicted  unit  plow  force  (Figure  4-3)  is  also  very  sensitive  to  the  surface  area  over  which  S22  is 
distributed.  This  bearing  area  is  slightly  truncated  near  the  tip  of  the  plow  in  the  present  plow/soil 
interaction  model  due  to  the  assumed  mode  of  failure  of  the  soil  at  the  tip.  Consequently,  the  unit 
resistance  force  may  tend  to  be  under-predicted. 

Time  histories  for  S3 3  are  presented  in  Figure  4-9e.  Steady-state  behavior  is  achieved  at 
all  centerline  elements  beyond  element  40.  The  gradual  rise  in  stress  level  while  the  plow  tip  is 
approaching  is  similar  to  S22;  both  S22  and  S33  sense  the  approaching  plow  tip,  where  Sll  does 
not.  Peak  values  of  the  S3  3  history  graphs  are  intermediate  between  Sll  and  S22  peak  values. 
Without  the  artificial  deformation  constraint  imposed  by  plane  strain  conditions,  S3 3  would 
likely  be  the  smallest  of  the  three  direct  stress  components  due  to  the  proximity  of  a  free  surface 
at  the  soil  bed.  In  Figure  4-9e,  large  gradients  of  S33  emanate  from  the  blade  surface  similar  to 
Sll  and  S22.  The  average  of  these  three  direct  stress  components  determines  the  strength  of  the 
soil  in  the  pressure-dependent  Drucker-Prager  soil  model.  If  S33  and  hence  the  soil  strength  are 
over-predicted,  the  unit  resistance  force  may  tend  to  be  over-predicted,  countering  the  under¬ 
predicting  effect  of  reduced  bearing  area. 
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Figure  4-9  Continued. 


4.4.2  Shear  Stress.  Shear  stresses  S12  are  shown  in  Figure  4-10.  Time  histories  rise  and 
fall  very  precipitously  at  the  centroids  of  centerline  elements  as  shown  in  Figure  4- 10a.  They  are 
nearly  flat  zero  due  to  the  imposed  symmetry  boundary  condition  along  the  centerline  until  the 
arrival  of  the  blade  tip.  A  slight  reversal  in  sign  is  apparent  before  and  after  arrival,  consistent 
with  the  reversal  noted  in  shear  strain  (Figure  4-8a).  Steady-state  behavior  is  evident  at  all 
centerline  elements  beyond  element  40.  Peak  values  in  shear  stresses  are  substantially  smaller 
than  peak  direct  stress  magnitudes.  Contours  of  shear  stress  are  shown  in  Figure  4- 10b.  Unlike 
direct  stress,  the  actual  maximum  shear  stress  tends  not  to  be  on  the  interface,  but  just  off  the 
surface  of  the  bevel  of  the  plow  tip.  Shear  stress  is  zero,  as  expected,  on  the  free  surface  formed 
aft  the  point  of  interface  separation. 

4.4.3  Mises  Equivalent  Stress.  Maximum  distortion  stress  is  perhaps  better  represented 
by  the  Mises  stress  shown  in  Figure  4-11,  an  invariant  measure  of  deviatoric  stress  behavior  in 
the  soil.  Stress  histories,  shown  in  Figure  4-1  la,  substantially  exceed  those  for  shear  stress 
component  S12  (Figure  4-10a).  Further,  they  sense  the  approaching  blade  with  a  gradual  increase 
in  value;  thereafter  a  precipitous  rise  and  fall  attends  actual  arrival  and  departure  of  the  plow  tip. 
Contour  data,  Figure  4-1  lb,  shows  the  maximum  MISES  stress  component  at  a  location  just  off 
the  blade  surface,  though  it  also  appears  to  be  related  to  the  shoulder  of  the  blade. 

4.4.4  Equivalent  Pressure  Stress.  The  equivalent  pressure  stress,  denoted  PRESS,  is 
presented  in  Figure  4-12.  PRESS  is  defined  as  the  negative  of  the  mean  effective  stress  in  these 
data,  so  it  can  be  easily  interpreted  as  the  latter  quantity  where  convenient.  PRESS  is  an  invariant 
stress  measure  related  to  the  dilatation  in  the  soil,  and  the  strength  of  the  soil  is  linearly  related  to 
PRESS  according  to  the  pressure-dependent  Drucker-Prager  soil  model.  Time  history  graphs  of 
PRESS  are  shown  in  Figure  4- 12a.  The  dilatant  material  expands  due  to  shearing  by  the  plow  as 
has  been  noted  in  the  discussion  of  soil  strain  (Figure  4-6).  The  mechanical  response  of  the 
partially  confined  soil  to  this  volume  expansion  is  a  state  of  compression  in  the  mean  effective 
stress,  as  the  plow  tip  arrives  and  departs.  The  pressure-dependent  strength  of  the  soil  at  a  point 
in  the  path  of  the  plow  along  the  centerline  correspondingly  rises  and  decays  with  the  arrival  and 
departure  of  the  plow  tip.  However,  when  observed  from  a  point  on  the  steadily  moving  plow  tip, 
the  plow  will  appear  continuously  confronted  with  soil  at  an  elevated  strength.  Figure  4-12b 
shows  steady-state  contour  values  of  PRESS.  They  can  be  interpreted  as  contours  proportional  to 
soil  strength,  and  show  soil  with  the  maximum  strength  bearing  on  the  surface  of  the  plow  tip, 
resisting  its  motion. 

Also,  while  the  time  histories  of  PRESS  and  MISES  are  apparently  very  similar  in  shape 
and  magnitude,  the  spatial  distributions  of  PRESS  and  MISES  differ  somewhat.  The  former  is 
similar  to  S 1 1 ,  S22,  and  S3  3,  all  of  which  are  maximum  on  the  surface  of  the  bevel  of  the  plow 
tip,  while  the  latter  is  maximum  just  off  this  surface. 
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(e)  Effective  stress  S3  3  histories  at  sampled  response  points  along  path  of  plow. 


(b)  Effective  stress  S12  distribution  near  plow  tip 
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(b)  Mises  effective  stress  distribution  near  plow  tip. 
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(b)  Mean  effective  stress  distribution  near  plow  tip. 


Soil  dilatancy  can  now  be  revisited.  If  element  volume  is  denoted  as  EVOL,  then  the 
change  in  element  volume  is  AEVOL  and  is  defined  as  the  algebraic  sum  of  the  direct  strains, 
El  1  and  E22.  Thus  AEVOL  has  an  elastic  and  a  plastic  component.  The  latter  component  is  the 
algebraic  sum  of  PEI  1  and  PE22,  and  the  former  component  is  the  mean  effective  stress  divided 
by  the  elastic  bulk  modulus  K  (5,000  psi  in  the  present  case).  Thus, 

AEVOL  =  El  1  +  E22 

=  (El  lel  +  E22el)  +  (El  IP1  +  E22P1) 

=  (-PRESS/K)  +  (PEI  1  +  PE22) 

The  volumetric  expansion  of  the  soil  can  be  determined  with  this  expression.  For  example,  where 
the  residual  value  of  PRESS  is  zero,  the  elastic  part  of  the  residual  volumetric  expansion  is  zero. 
The  residual  expansion  of  the  material  in  the  wake  is  therefore  the  sum  of  residual  plastic  strain 
components  PEI  1  and  PE22  (Figures  4-6c  and  h).  In  this  case,  the  maximum  residual  volumetric 
expansion  equals  27  percent,  as  before. 

4.4.5  Summary  of  Model  Stresses  in  the  Soil.  The  mechanical  response  of  the  soil 
skeleton  is  characterized  by  effective  stresses.  Like  the  strain  response,  the  effective  stress 
response  is  primarily  confined  to  one  or  two  blade  widths  from  the  plow  tip,  and  stress  gradients 
are  very  large.  Effective  stress  components  rise  and  fall  precipitously  during  arrival  and 
departure  of  the  plow  tip,  as  do  their  corresponding  strain  components.  However,  unlike  strains, 
stresses  return  to  zero  as  the  soil  is  unloaded  and  there  is  no  residual  stress  response.  Although 
plow  tip  forces  are  based  on  total  stress,  the  distribution  of  effective  stress  over  the  plow  tip  is 
consistent  with  these  forces.  The  assumed  failure  mechanism  for  cleaving  the  soil  at  the  plow  tip 
causes  a  small  reduction  in  the  area  of  contact  at  the  tip  of  the  plow  which  is  countered  by  an 
over-prediction  of  soil  strength.  Further  features  of  the  computed  stress  response  in  the  plow/soil 
interaction  model  are: 

•  The  location  of  the  maximum  distortion  stress  corresponds  to  the  shoulder  of  the  plow 
tip. 

•  The  distribution  of  the  lateral  stress  component  over  the  bevel  of  the  plow  tip  is 
approximately  linear,  varying  from  zero  at  the  tip  to  a  maximum  at  the  shoulder.  This 
result  is  consistent  with  wear  patterns  observed  on  steel  blades  employed 
continuously  in  the  field  which  show  maximum  wear  on  the  shoulder  of  the  plow  tip. 

•  The  distribution  of  the  in-line  component  of  stress  over  the  bevel  of  the  plow  tip  is 
more  uniform,  although  slightly  truncated  at  the  tip. 

•  Soil  strength  is  proportional  to  mean  pressure  stress  according  to  the  Drucker-Prager 
soil  model  employed.  Contours  of  pressure  stress  demonstrate  how  the  plow  tip,  to 
progress  forward,  must  continuously  overcome  soil  at  maximum  strength  conditions. 
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4.5  Pore  Pressure  in  the  Soil 


Whereas  effective  stress  reflects  the  behavior  of  the  soil  skeleton  or  solid  phase  of  the 
soil,  pore  pressure  reflects  the  behavior  of  the  fluid  phase  of  the  soil.  It  is  denoted  by  the  variable 
POR.  Unlike  effective  stress,  pore  pressure  is  computed  and  reported  at  the  node  points;  no 
averaging  or  interpolation  occurs  in  post  processing  of  these  data. 

In  Figure  4-13,  a  negative  or  reduced  pore  pressure  field  is  shown  which  results  when  a 
dilatant  soil  is  subjected  to  the  shearing  action  of  the  plow.  This  response  evidences  the  coupling 
that  exists  between  the  fluid  and  solid  phases  of  the  soil.  The  maximum  negative  pore  pressure 
coincides  with  the  plow  tip  where  maximum  volumetric  expansion  occurs. 

The  pore  pressure  field  is  more  diffuse  in  time  and  space  than  the  effective  stress  field. 
While  the  absence  of  large  gradients  facilitates  accurate  modeling  of  the  pore  pressure  field  with 
low  order  finite  elements,  its  diffuse  nature  makes  it  more  difficult  to  simulate  without 
contaminating  the  solution  in  the  region  of  interest  near  the  plow  tip  due  to  effects  from  the 
artificial,  finite  boundaries  of  the  model.  Nonetheless,  there  are  several  similarities  in  the 
responses  of  computed  pore  pressure  and  measured  pore  pressure  (with  the  laboratory  apparatus. 
Figure  1-1). 

In  Figure  4-1 3a,  time  histories  of  pore  pressures  are  presented  at  seven  uniformly  spaced 
nodes  along  the  centerline  (Figure  4-2).  Pore  pressure  change  is  clearly  sensed  well  in  advance  of 
the  plow,  a  response  also  observed  experimentally.  There  is  a  substantial  overshoot  in  peak  pore 
pressure  as  the  plow  begins  to  penetrate  into  the  soil  during  the  startup  transient  phase  of  the 
motion.  According  to  these  data,  pore  pressure  response  becomes  steady-state  when  the  plow  has 
penetrated  to  node  101  and  beyond  in  the  plow/soil  interaction  model.  The  peak  steady-state  pore 
pressures  shown  are  representative  of  measurements  taken  near  (but  not  on)  the  centerline  in  the 
laboratory  apparatus  for  the  particular  plow  configuration  and  plow  velocity  corresponding  to 
this  simulation.  In  Figure  4-1 3b,  the  maximum  negative  pore  pressure  or  trough,  envelops  the  tip 
of  the  plow  blade.  The  trough  is  roughly  circular  in  shape  and  surrounds  the  leading  edge  of  the 
blade,  both  of  which  are  observations  made  for  measured  pore  pressure  distribution  in  the 
laboratory  apparatus.  Also,  refer  to  Rutter  and  V oss  ( 1 994)  where  similar  pore  pressure  contours 
were  reconstructed  from  measured  data. 

The  reduced  view  in  Figure  4-1 3c  shows  the  diffuse  nature  of  pore  pressure  over  the 
entire  field  of  the  plow/soil  interaction  model.  The  graph  of  pore  pressure  distribution  along  the 
path  of  the  plow  shows  the  peak  response  at  the  plow  tip,  and  the  gradual  decay  forward  and  aft 
of  the  plow  tip.  Pore  pressure  gradients  are  relatively  small  in  comparison  to  gradients  of 
effective  stress  and  strain  in  the  soil  skeleton.  However,  where  the  pore  pressure  does  not 
approach  zero  at  a  boundary  of  the  plow/soil  interaction  model  (except  at  the  centerline),  the  pore 
fluid  region  may  be  considered  artificially  truncated,  and  the  response  may  not  be  an  accurate 
representation  of  pore  pressure  in  that  region. 
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(a)  Pore  pressure  histories  at  sampled  response  points  along  path  of  plow 
Figure  4-13.  Pore  pressure  in  soil  at  5  fps  steady-state  plow  velocity. 
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(b)  Pore  pressure  distribution  near  plow  tip. 
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(c)  Pore  pressure  field  for  plow/soil  interaction  model. 
Figure  4-1 3  Continued. 


In  summary,  the  pore  pressure  field  computed  by  the  plow/soil  interaction  model  is  very 
diffuse  as  contrasted  with  the  effective  stress  field  computed  by  the  model  which  is  characterized 
by  large  gradients. 

•  The  pore  pressure  field  computed  by  the  plow/soil  interaction  model  is  everywhere 
negative,  behavior  which  is  expected  of  a  dilatant  soil  subjected  to  shear  stress. 

•  Pore  pressure  magnitudes  are  representative  of  measured  pore  pressures  in  the 
laboratory  apparatus.  In  comparison  to  maximum  effective  stresses,  the  maximum 
pore  pressures  are  small,  explaining  why  at  high  stress  locations  the  total  and 
effective  stress  responses  appear  equal. 

•  The  roughly  circular  trough  in  the  contour  map  of  computed  pore  pressure  envelops 
the  plow  tip.  It  corresponds  to  the  area  of  maximum  volumetric  expansion,  and  the 
location  of  maximum  negative  pore  pressure  measured  in  the  laboratory  apparatus. 
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5.0  PARAMETER  STUDY 


This  parameter  study  demonstrates  the  kind  of  engineering  data  that  can  be  derived  from 
a  rational  numerical  model  of  plow/soil  interaction.  The  parameters  chosen  relate  to  properties  of 
the  plow;  they  are  the  plow  velocity,  the  coefficient  of  friction  between  the  blade  and  soil,  and 
the  shape  of  the  blade  as  represented  by  the  included  angle  of  the  blade  tip.  Baseline  conditions 
or  standard  blade  properties  were  selected  as  follows:  10  fps,  no  blade  or  interface  friction,  and  a 
1 /2-inch-thick  blade  with  a  60-degree  blade  tip  included  angle  (Figure  1-lb).  The  effect  that 
velocity,  friction  and  shape  parameters  have  on  plow  forces,  and  on  soil  stress  and  pore  pressure 
along  the  path  of  travel  are  studied. 

The  effect  of  velocity  on  tow  force  can  be  studied  more  accurately  by  experiment  in  the 
portable  test  apparatus  (Figure  1-la)  since  the  present  two-dimensional  numerical  model  predicts 
only  unit  tow  force,  not  total  force.  Therefore,  no  attempt  is  made  to  directly  compare 
experimental  and  numerical  data,  although  frequently  experimental  observations  and  model 
predictions  are  qualitatively  compared  in  the  discussion.  Conversely,  the  effect  of  parameters  like 
blade  friction  and  blade  shape  can  perhaps  be  studied  more  easily  with  the  numerical  model. 
These  parameters  have  not  been  studied  experimentally  in  the  portable  test  apparatus  to  date.  The 
study  demonstrates  how  the  numerical  model  can  complement  experimental  approaches  to  gain 
understanding  of  plow/soil  interaction. 

5.1  Effect  of  Plow  Velocity  on  Performance  and  Soil  Response 

The  steady-state  velocities  listed  in  Table  3-2  are  used  to  study  the  effect  of  plow 
velocity.  The  plow/soil  interaction  model  was  run  for  each  of  these  velocities  until  the  plow 
achieved  steady-state  conditions  and  traveled  approximately  the  same  embedment  distance 
(Figure  4-1).  A  frictionless  blade  was  assumed  and  the  standard  1 /2-inch-thick  blade  with  a  60- 
degree  blade  tip  angle  was  assumed. 

5.1.1  Plow  Forces.  Histories  of  tow  and  lateral  force  components  acting  on  the  plow  are 
shown  in  Figure  5-1  for  the  case  of  10  fps  steady-state  velocity.  These  are  forces  per  unit  depth 
of  the  plow  blade,  or  unit  plow  force  components.  The  steady-state  forces  for  each  component 
are  obtained  graphically  by  averaging  the  steady  portion  of  the  history  data  shown.  The  lower 
curve  is  the  tow  force  component  for  half  the  blade  thickness,  so  the  tow  force  on  the  full  1/2- 
inch-thick  blade  is  twice  the  steady-state  value  shown  in  the  graph.  This  clearly  makes  it 
approximately  equal  to  the  lateral  component  shown  in  the  upper  graph. 

The  linear  relationship  that  these  plow  force  components  apparently  have  with  steady- 
state  plow  velocity,  over  the  range  of  velocities  studied,  is  shown  in  Figure  5-2.  This  is  in 
agreement  with  experimental  observation.  However,  there  is  little  evidence  of  a  critical  velocity 
in  these  data  such  as  has  been  observed  experimentally  (Figure  1-2),  and  as  proposed  by  the 
semi-empirical  model  (Figure  2-1).  This  is  perhaps  related  to  the  numerical  model’s  inability  to 
allow  softening  in  the  Drucker-Prager  model  of  the  solid  phase,  or  to  allow  a  cavitation  pressure 
cutoff  in  modeling  the  fluid  phase.  On  closer  inspection,  a  quadratic  least  squares  curve  fit  to  the 
tow  force  data  is  perceptively  better  than  a  linear  curve  fit  indicating  only  a  very  slight  tendency 
toward  critical  velocity  behavior. 
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Figure  5-1 .  Unit  plow  force  histories  at  10  fps  steady-state  plow  velocity. 


Unit  Plow  Force  (lb/ in) 


A  static  tow  force,  i.e.,  the  tow  force  at  zero  velocity,  can  be  estimated  from  this  graph 
using  the  method  of  Section  4.2.3.  The  static  tow  force  is  typically  non-zero,  as  shown  in  the 
experimental  data  (Figure  1-lb).  However,  comparison  between  computed  unit  plow  forces  and 
measured  total  plow  forces  can  be  misleading  for  the  two-dimensional  plane  strain  model  cannot 
predict  total  force. 

The  apparent  equality  between  the  unit  tow  and  lateral  force  components  occurs  for  the 
entire  range  of  velocities  studied  for  the  case  of  the  standard  ffictionless  blade.  There  is  no 
apparent  reason  why  this  should  be  so,  but  it  is  a  useful  rule  of  thumb  since  lateral  force  has  not 
been  measured. 

Unit  plow  force  is  an  integral  of  contact  force  distribution  over  the  tip  of  the  plow  blade, 
as  shown  for  the  10  fjps  case  in  Figure  5-3.  Contact  forces  are  forces  per  unit  length  (lb/in.)  along 
the  plow/soil  interface.  Typically,  the  predicted  profile  of  the  lateral  contact  force  along  the 
interface  varies  nearly  linearly  from  zero  forward  of  the  tip  of  the  blade,  to  a  maximum  at  the 
shoulder  of  the  blade,  and  then  rapidly  to  zero  just  aft  of  the  shoulder  where  interface  separation 
begins.  The  maximum  contact  force  would  be  expected  to  increase  linearly  with  plow  velocity. 

5.1.2  Behavior  of  Principal  Stresses  in  the  Soil.  Total  stresses  in  the  soil  which  act  at 
the  plow/soil  interface  account  for  the  aforementioned  contact  forces  and  the  overall  resistance 
the  plow  experiences  during  penetration.  Total  stress  is  the  algebraic  sum  of  the  effective  stress 
and  pore  pressure  in  the  soil,  but  the  former  is  generally  much  larger  than  the  latter.  A  reasonably 
complete  picture  of  resistance  is  therefore  provided  by  the  distribution  of  principal  (effective) 
stresses  near  the  plow  tip. 

Contours  of  maximum  and  minimum  principal  stresses  in  the  soil  near  the  tip  of  the  plow 
are  presented  in  Figure  5-4  for  the  case  of  the  standard  frictionless  plow  blade  moving  at  10  fps 
(intermediate  principal  stress  is  omitted).  The  maximum  shear  stress  is  also  indirectly  included  in 
these  data,  and  is  obtained  as  one-half  the  difference  between  the  minimum  (SP1,  lowest 
negative  or  compression  stress)  and  maximum  (SP3,  highest  negative  or  compression  stress) 
principal  stresses  presented  in  Figures  5-4a  and  b,  respectively.  The  greatest  stress  acts  directly 
on  the  shoulder  of  the  blade  where  field  observation  correspondingly  reveals  most  of  the  wear. 

Variation  of  principal  stresses  with  steady-state  velocity  is  illustrated  in  Figure  5-5.  Data 
plotted  are  the  extreme  contour  values  (Figure  5-4)  for  each  velocity.  A  least  squares  2nd  order 
polynomial  curve  fit  to  these  data  is  clearly  superior  to  a  linear  curve  fit.  Principal  stresses 
therefore  appear  to  be  approaching  a  lower  bound  or  saturation  level.  Also,  the  minimum  and 
maximum  curves  draw  closer  (by  about  37  percent  in  this  case)  over  the  range  of  velocities 
studied.  Half  their  difference,  at  a  particular,  velocity,  is  an  upper  bound  on  the  maximum  shear 
stress  that  exists  near  the  plow  tip  at  that  velocity.  Accordingly,  the  maximum  shear  stress 
decreases,  and  the  two  principal  stress  components  increase  (negatively)  toward  a  bound,  with 
increasing  steady-state  plow  velocity. 

5.1.3  Pore  Pressure  Behavior.  Due  to  the  shearing  action  of  the  plow  tip,  the  dilatant 
soil  expands  as  the  plow  approaches,  and  because  of  the  plasticity  in  this  deformation,  element 
volumes  in  the  path  of  the  plow  (along  the  centerline,  CL,  of  the  model)  increase  permanently  to 
the  residual  levels  shown  in  Figure  5-6,  for  the  standard  frictionless  blade  moving  at  10  fps.  Pore 
water  migration  into  the  void  is  thereby  initiated  according  to  soil  permeability,  a  prescribed 
function  of  void  ratio  (Figure  3-12).  The  corresponding  pore  pressure  histories  are  shown  in 
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Figure  5-la.  Data  sampled  at  node  points  above  node  61  (Figure  4-2),  vary  similarly  and  verify 
that  steady-state  behavior  exists  in  the  fluid  phase  of  the  plow/soil  interaction  model.  Pore 
pressures  are  everywhere  negative,  and  the  minimum  steady-state  pore  pressure  can  be 
determined  from  these  data.  The  minimum  pressure  or  trough  in  pore  pressure  at  the  plow  tip 
created  by  shearing  of  the  dilatant  soil  is  evident  in  Figure  5-7b.  The  algebraic  sum  of  this 
negative  pressure  and  the  effective  stress  (compressive,  Figure  5-4a)  near  the  tip,  i.e.,  the  total 
stress,  constitutes  the  source  of  resistance  to  plow  penetration.  According  to  these  results,  the 
minimum  pore  pressure  engulfing  the  blade  tip  an  order  of  magnitude  smaller  than  the  principal 
stress  at  the  shoulder  of  the  blade  tip. 

The  variation  of  minimum  pore  pressure  at  the  blade  tip  with  steady-state  velocity  is 
shown  to  be  linear  in  Figure  5-8.  The  minimum  pore  pressure  decreases  without  bound  because 
the  plow/soil  interaction  model  does  not  include  any  provision  for  cavitation  in  the  pore  fluid. 

5.2  Effect  of  Plow  Friction  on  Performance  and  Soil  Response 

The  smoothness/roughness  condition  of  the  plow  blade  is  undoubtedly  an  important 
parameter  as  it  relates  to  performance,  effectiveness,  and  operational  cost  of  trenching 
operations.  It  can  be  represented  by  a  kinetic  coefficient  of  friction  at  the  plow/soil  interface  in 
conjunction  with  the  standard  Coulomb  friction  model  in  the  plow/soil  interaction  model.  The 
aforementioned  parametric  results  (Figures  5-1  through  5-8)  pertain  to  a  frictionless  interface 
between  the  plow  and  soil.  In  this  parameter  study  the  coefficient  of  friction  between  the  blade 
and  soil  was  varied  to  determine  its  effect  on  plow  forces  and  on  behavior  of  the  soil.  The  effect 
of  small  amounts  of  friction,  i.e.,  5  and  10  percent,  are  studied  with  the  standard  blade  moving  at 
10  fps  steady-state  velocity.  The  perfectly  smooth  case  (zero  percent)  is  considered  the  baseline 
case. 


5.2.1  Plow  Forces.  Unit  plow  force  histories  including  the  effect  of  friction  are 
presented  in  Figure  5-9.  Predictably,  both  components  of  plow  force  have  increased  compared  to 
the  baseline  case  (Figure  5-1),  and  both  components  increase  with  the  value  of  the  coefficient  of 
friction  as  seen  by  comparing  the  graphs  of  Figures  5-9a  and  b  which  show  the  results  for  5  and 
10  percent  friction,  respectively.  When  friction  is  present,  the  lateral  force  is  somewhat  less  than 
the  tow  force,  and  the  aforementioned  equality  between  lateral  and  tow  force  components 
apparently  pertains  only  to  the  case  of  a  frictionless  interface  (i.e.,  the  lower  curves  have  drawn 
closer  to  the  upper  curves  as  friction  is  increased,  and  when  multiplied  by  two  to  correct  for  a 
half-blade  thickness,  they  exceed  the  magnitudes  of  the  upper  curves). 

This  is  seen  more  clearly  in  Figure  5-10,  where  the  tow  and  lateral  force  components  are 
shown  to  vary  with  the  coefficient  of  friction  due  to  the  Coulomb  friction  sliding  interface 
model.  The  normalized  values  shown  are  obtained  by  dividing  the  two  steady-state  plow  force 
components  by  their  respective  baseline  values,  which  were  essentially  equal.  The  tow  force 
component  increases  linearly  and  more  rapidly  with  blade  friction  than  the  lateral  force 
component.  Indeed,  the  plow/soil  interaction  model  predicts  that  the  tow  force  increases  35 
percent  for  a  coefficient  of  friction  of  10  percent  between  the  blade  and  soil,  for  the  case  of  the 
standard  blade  moving  with  10  fps  steady-state  velocity.  The  lateral  force  increases  less  than 
linearly  towards  an  apparent  bound  of  only  4  percent. 
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Figure  5-3.  Distribution  of  lateral  interface  contact  force  on  plow  at  10  fps 
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(a)  Minimum  principal  stress  SP1  distribution  near  plow  tip. 

Figure  5-4.  Principal  stress  components  in  soil  at  10  fps  steady-state  plow  velocity. 
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(b)  Maximum  principal  stress  SP3  distribution  near  plow  tip. 
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Figure  5-5.  Effect  of  plow  velocity  on  principal  stresses  in  soil. 
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Figure  5-6.  Element  volume  change  histories  at  10  fps  steady-state  plow  velocity. 
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(a)  Pore  pressure  histories  at  sampled  response  points  along  path  of  plow. 
Figure  5-7.  Pore  pressure  in  soil  for  10  fps  steady-state  plow  velocity. 
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Figure  5-8.  Effect  of  plow  velocity  on  minimum  pore  pressures. 
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The  corresponding  steady-state  lateral  contact  forces  on  the  tip  of  the  blade  vary  along  the 
interface  as  shown  in  Figure  5-11,  i.e.,  they  vary  nearly  linearly  from  zero  at  the  tip  to  a 
maximum  at  the  shoulder  of  the  blade,  as  they  did  for  the  baseline  case  (Figure  5-3).  The 
maximum  value  of  the  lateral  contact  force  for  10  percent  friction  is  1 1.5  percent  greater  than  for 
the  case  of  no  friction. 

5.2.2  Behavior  of  Principal  Stresses  in  the  Soil.  The  two  steady-state  principal  stress 
components  for  the  case  of  5  and  10  percent  friction  are  plotted  in  Figures  5-12a  and  b,  and  5-12c 
and  d,  respectively.  Contours  of  their  spatial  distribution  focus  on  the  shoulder  of  the  plow  blade. 
When  compared  to  the  baseline  case  (Figure  5-4)  it  is  seen  that  the  effect  of  friction  is  to 
decrease  the  minimum  principal  stress  and  increase  the  maximum  principal  stress.  For  example, 
with  the  standard  blade  moving  at  10  fps,  adding  10  percent  friction  decreases  the  lowest 
minimum  stress  by  28  percent,  and  increases  the  highest  maximum  stress  by  22  percent.  The  two 
extreme  principal  stresses,  normalized  on  their  average  for  the  baseline  case  (-735  psi),  are 
graphed  against  the  coefficient  of  friction  in  Figure  5-13.  Half  the  difference  in  the  two  extreme 
curves  shown  approximates  the  maximum  shear  stress  near  the  plow  tip.  So  the  maximum  shear 
stress  increases  dramatically  with  increasing  friction.  Indeed,  by  adding  1 0  percent  friction  in  the 
plow/soil  interaction  model  the  maximum  shear  stress  increases  by  a  factor  of  3.0.  Thus,  a 
rougher  blade  surface  develops  substantially  more  shear  stress  in  the  soil. 

5.2.3  Pore  Pressure  Behavior.  Histories  of  volumetric  expansion  of  the  soil  in  the  path 
of  the  plow  for  the  cases  of  5  and  10  percent  friction  are  shown  in  Figures  5- 14a  and  b, 
respectively.  Only  a  small  increase,  due  to  blade  friction,  is  noted  in  comparison  to  the  baseline 
case  (Figure  5-6).  By  contrast,  it  will  be  shown  in  the  shape  parameter  study  discussed  below, 
that  a  dull  blade  tip  causes  a  substantial  increase  in  volumetric  expansion  of  the  dilatant  soil. 

Pore  pressure  histories  and  contours  of  pore  pressure  distribution  for  the  cases  of  5  and  1 0 
percent  friction  are  shown  in  Figures  5- 15a  and  b,  and  Figures  5- 15c  and  d,  respectively. 
Correspondingly,  when  compared  to  the  baseline  case  (Figure  5-7),  there  is  no  noticeable 
difference  due  to  blade  friction.  Only  on  very  close  inspection  of  the  minimum  steady-state 
values  in  the  pore  pressure  history  data,  can  a  very  slight  effect  be  seen.  This  is  shown  in  Figure 
5-16,  where  minimum  pore  pressures,  normalized  by  the  corresponding  value  for  the  baseline 
case,  are  graphed  against  the  coefficient  of  friction.  Adding  10  percent  friction  lowers  the 
minimum  pore  pressure  by  only  4  percent.  This  is  consistent  with  the  small  effect  noted  in 
volumetric  expansion. 

5.3  Effect  of  Plow  Tip  Shape  on  Performance  and  Soil  Response 

The  effectiveness  of  the  blade  in  cutting  through  soil  is  important,  and  is  related  to  the 
profile  or  shape  of  the  blade  tip.  The  shape  of  the  blade  tip  has  an  influence  on  the  towing  force 
of  blades  as  indicated  with  the  single-step  finite  element  model  discussed  earlier  (Figure  2. 2d). 
Further,  the  shape  of  the  blade  tip  may  be  altered  due  to  wear  on  blades  used  continuously  in  the 
field,  causing  a  change  in  their  performance.  For  example,  it  has  been  observed  that  the  shoulder 
of  the  blade  tip  wears,  so  that  lower  sections  of  the  blade  entry  are  visibly  sharper,  although  the 
very  tip  of  the  lower  blade  itself  is  somewhat  more  dulled.  This  parameter  study  helps  to 
understand  the  corresponding  consequence  of  continued  blade  usage. 
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(a)  Five  percent  friction. 

Figure  5-1 1 .  Distribution  of  lateral  interface  contact  stress  on  plow  at  10  fps  steady-state  plow  velocity 
with  friction. 
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(b)  Ten  percent  friction. 
Figure  5-11  Continued. 


(N  Cl  IN  M 


N  (N  fN  fN 

o  o  -  - 


o  o  © 


in  w  a  « 

t-i  in  «•  (M  « 


W  N  W  W  W 


rt  M  M  H 


118 


(a)  Minimum  principal  stress  SP1  distribution  near  plow  tip  with  5  percent  friction. 

Figure  5-12.  Principal  stress  components  in  soil  at  10  fps  steady-state  plow  velocity  with 
blade  friction. 
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(b)  Maximum  principal  stress  SP3  distribution  near  plow  tip  with  5  percent  friction. 

Figure  5-12  Continued. 
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(d)  Maximum  principal  stress  SP3  distribution  near  plow  tip  with  10  percent  friction. 

Figure  5-12  Continued. 
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(a)  Five  percent  friction. 

Figure  5-14.  Element  volume  change  histories  at  10  fps  steady-state  plow  velocity  with  blade 
friction. 


(b)  Ten  percent  friction. 


(d)  Pore  pressure  distribution  in  soil  surrounding  plow  with  10  percent  friction. 

F igure  5-15  Continued. 
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Figure  5-16.  Effect  of  blade  friction  on  minimum  pore  pressure. 


In  the  parameter  study  the  effect  of  the  geometry  of  the  blade  tip  is  studied  by  exercising 
the  plow/soil  interaction  model  for  three  cases;  a  sharp  tip,  a  conventional  or  standard  tip,  and  a 
dull  tip,  defined  by  included  tip  angles  of  30,  60,  and  90  degrees,  respectively.  The  standard 
frictionless  blade  with  the  60  degree  included  tip  angle  (Figure  1-lb),  moving  at  10  fps  is  the 
baseline  case. 

5.3.1  Plow  Forces.  Unit  plow  force  histories  for  sharp  and  dull  frictionless  blades 
moving  at  10  fps  are  shown  in  Figure  5-17.  Steady-state  plow  forces  are  clearly  very  different 
from  those  for  the  baseline  case  (Figure  5-1),  where  the  tow  and  lateral  force  components  were  in 
fact  nearly  equal.  The  tow  force  for  the  sharp  blade,  shown  in  Figure  5- 17a,  is  significantly 
lower,  while  the  lateral  force  is  only  somewhat  lower  than  in  the  baseline  case.  The  lateral  force 
is  about  twice  the  tow  force  for  the  sharp  blade  case.  The  tow  force  for  the  dull  blade,  shown  in 
Figure  5-1 7b,  is  significantly  higher,  while  the  lateral  force  is  only  somewhat  higher  than  the 
baseline  case.  The  lateral  force  is  substantially  less  than  the  tow  force  for  the  dull  blade  case. 

These  data  are  summarized  in  Figure  5-18,  where  the  tow  and  lateral  forces,  normalized 
on  their  (nearly  equal)  baseline  values,  are  graphed  against  the  blade  tip  included  angle.  The 
lateral  force  component  is  insensitive  but  the  tow  force  component  is  very  sensitive  to  blade  tip 
angle.  For  the  sharp  blade,  the  tow  force  is  reduced  to  approximately  45  percent  of  that  for  the 
standard  blade.  For  the  dull  blade,  the  tow  force  is  increased  to  approximately  70  percent  of  that 
for  the  standard  blade.  The  blade  tip  included  angle  has  a  substantial  influence  on  towing  force 
while  affecting  the  lateral  force  very  little. 

The  distribution  of  lateral  contact  forces  along  the  plow/soil  sliding  interface  for  the  two- 
blade  configurations  is  shown  in  Figure  5-19.  As  expected,  contact  forces,  i.e.,  forces  per  unit 
length  along  the  contact  surface,  vary  nearly  linearly  from  zero  at  the  tip  to  a  maximum  at  the 
shoulder  of  the  blade.  While  the  maximum  contact  forces  for  the  two-blade  configurations  shown 
and  for  the  baseline  case  (Figure  5-3)  differ  substantially,  the  areas  beneath  all  three 
distributions,  i.e.,  the  total  lateral  contact  forces,  differ  only  modestly.  Assuming  a  linear 
distribution,  the  areas  under  the  three  contact  force  distributions  shown  all  agree  within  7  percent 
with  the  steady-state  lateral  plow  forces  obtained  from  force  history  data.  This  is  consistent  with 
the  above  finding  that  the  lateral  component  of  the  plow  force  is  insensitive  to  blade  tip  angle 
(Figure  5-18).  However,  because  the  maximum  lateral  contact  force  increases  substantially  with 
blade  tip  angle  (Figure  5-19),  sharper  blade  tips  would  not  sustain  as  much  shoulder  wear  as 
duller  blade  tips. 

5.3.2  Behavior  of  Principal  Stresses  in  the  Soil.  Principal  stress  distributions  in  the 
soil  for  the  sharp  and  dull  blades  are  shown  in  Figure  5-20.  In  the  neighborhood  of  the  plow  tip 
these  effective  stresses  are  compressive,  and  they  are  focused  on  the  shoulder  of  the  tip  of  the 
blade.  As  expected,  sharper  blades  induce  less  stress  in  the  soil  than  duller  blades.  These  stresses, 
compared  with  those  for  the  standard  blade  (Figure  5-4),  shows  that  the  minimum  principal  stress 
increases  (becomes  less  compressive)  by  48  percent  for  the  sharp  blade,  and  decreases  (becomes 
more  compressive)  64  percent  for  the  dull  blade.  The  variation  of  minimum  principal  stress  is 
nearly  linear  with  blade  tip  angle.  The  maximum  principal  stress  increases  60  percent  for  the 
sharp  blade,  but  only  decreases  8  percent  for  the  dull  blade.  The  variation  of  maximum  principal 
stress  saturates,  i.e.,  approaches  a  horizontal  asymptote.  This  variation  can  be  seen  in  Figure  5-21 
where  the  two  principal  stress  components  are  normalized  on  the  average  principal  stress  for  the 
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standard  blade  (-735  psi),  and  graphed  against  blade  tip  angle.  Saturation  of  maximum  principal 
stress  beyond  an  angle  of  60  degrees  causes  a  large  increase  in  maximum  shear  stress  (difference 
between  curves)  for  the  dull  blade.  This  is  compared  with  a  modest  23  percent  decrease  in 
maximum  shear  stress  for  the  sharp  blade.  At  the  same  velocity,  the  dull  blade  tip  would  be  far 
more  effective  in  overcoming  the  shear  strength  of  the  soil  and  producing  a  trench  but  at  the  cost 
of  a  higher  tow  force  (Figure  5-18). 

5.3.3  Pore  Pressure  Behavior.  The  shape  of  the  blade  tip  has  a  pronounced  effect  on 
volumetric  expansion  of  the  dilatant  soil  along  the  path  of  the  plow  as  shown  in  Figure  5-22. 
Comparing  the  residual  volumetric  expansions  for  the  sharp  and  dull  blades  in  Figures  5-22a  and 
b,  respectively,  with  the  residual  volumetric  expansion  for  the  baseline  case  (Figure  5-6),  shows 
noticeably  less  expansion  occurring  with  the  sharp  blade  and  noticeably  more  expansion 
occurring  for  the  dull  blade.  Thus,  while  friction  may  have  had  a  small  influence,  blade  tip  shape 
has  a  significant  influence  on  residual  volumetric  expansion  of  soil  in  the  path  of  the  plow. 

Correspondingly,  the  pore  pressure  histories  and  pore  pressure  contours  would  be 
expected  to  vary  significantly  with  the  shape  of  the  blade,  and  these  data  are  shown  in  Figure 
5-23.  The  significant  pore  pressure  is  the  minimum  steady-state  value  shown  in  the  history  data, 
and  also  shown  as  a  pressure  trough  engulfing  the  blade  tip  in  the  contour  data  (history  data  are 
somewhat  more  accurate).  The  variation  of  minimum  pore  pressure  with  blade  tip  angle  is  shown 
in  Figure  5-24  where  pressures  have  been  normalized  on  the  minimum  value  for  the  baseline  case 
(-17  psi).  Thus  a  30  degree  blade  tip  angle  will  increase  the  minimum  pore  pressure  by 
approximately  30  percent,  and  a  90  degree  blade  tip  angle  will  decrease  the  minimum  pore 
pressure  approximately  1 8  percent. 
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(a)  Lateral  and  tow  (half)  force  components  per  unit  depth  of  blade  with  sharp  tip. 

Figure  5-17.  Unit  plow  force  histories  at  10  fps  steady  state  plow  velocity  with  sharp  and  dull 
blades. 
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(b)  Lateral  and  tow  (half)  force  components  per  unit  depth  of  blade  with  dull  tip 

Figure  5-17  Continued. 
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Figure  5-18.  Effect  of  blade  tip  angle  on  steady-state  plow  forces. 
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Figure  5-19.  Distribution  of  lateral  interface  contact  forces  on  plow  at  10  fps 
steady-state  plow  velocity  with  sharp  and  dull  blades. 
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Figure  5-19  Continued 


(a)  Minimum  principal  stress  SP1  distribution  near  sharp  tip. 

Figure  5-20.  Principal  stress  components  in  soil  at  10  fps  steady-state  plow  velocity  with 
sharp  and  dull  blades. 
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(b)  Maximum  principal  stress  SP3  distribution  near  sharp  tip. 
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(c)  Minimum  principal  stress  SP1  distribution  near  dull  tip. 
Figure  5-20  Continued. 
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(d)  Maximum  principal  stress  SP3  distribution  near  dull  tip. 
Figure  5-20  Continued. 


20  30  40  50  60  70  80  90  100 


Blade  Tip  Included  Angle  (deg.) 


Figure  5-21 .  Effect  of  blade  tip  angle  on  principal  stresses  in  soil. 
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Figure  5-22  Continued. 
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(b)  Pore  pressure  distribution  in  soil  surrounding  plow  with  sharp  blade 

Figure  5-23  Continued. 
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(c)  Pore  pressure  histories  at  sampled  response  points  along  path  of  plow  with  dull  blade. 
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Figure  5-24.  Effect  of  blade  tip  angle  on  minimum  pore  pressure. 


6.0  SUMMARY  AND  CONCLUSIONS 


The  development  of  a  rational  numerical  model  for  the  interaction  of  a  thin  blade  plow 
with  seafloor  sand  during  steady-state  trenching  in  support  of  cable  burial  operations  was 
addressed.  Earlier  attempts  to  explain  differences  between  measured  tow  forces  in  the  field  and 
predicted  tow  forces  based  on  conventional  theories  were  unsuccessful.  More  recent  laboratory 
experiments  at  NFESC  have  overcome  difficulties  in  preparation  of  sands  that  replicate  100 
percent  saturation  conditions  in  seafloor  soil,  and  have  successfully  measured  tow  forces  on  scale 
models  of  thin  blade  plows  moving  with  constant  velocity.  The  data  showed  that  the  relationship 
between  tow  force  and  velocity  was  linear  below  a  critical  velocity,  beyond  which  it  leveled  off 
to  a  constant  tow  force.  Development  of  numerical  models  was  also  authorized  to  investigate  the 
mechanics  of  plow/soil  interaction,  to  compute  plow  forces,  to  provide  a  rational  analysis  tool  for 
analyzing  the  effects  of  alternative  seafloor  soils  and  plow  design  parameters,  and  to  assess  the 
state-of-the-art  in  commercial  finite  element  technology  relative  to  these  objectives. 

Numerical  modeling  of  plow/soil  interaction  is  a  complex  problem  irrespective  of 
whether  state-of-the-art  commercial  finite  element  technology  or  existing  advanced  special 
purpose  software  is  employed  in  the  solution.  The  fundamental  problem  is  to  model  the  coupled 
governing  field  equations  for  saturated  soil  in  the  presence  of  a  rigid  body  having  sliding 
interfaces  while  it  is  penetrating  seafloor  soil.  The  field  equations  are  those  for  equilibrium  of  an 
elasto-plastic  porous  medium,  coupled  to  those  for  pore  water  flow  within  the  medium.  Effective 
stress  analysis  is  required,  friction  on  sliding  interfaces  must  be  provided  for,  and  the  stress  state 
in  the  soil  is  such  that  the  analysis  ideally  should  be  carried  out  in  three  dimensions. 

The  state  of  commercial  software  coerces  a  multistep  transient  solution  approach  for  the 
steady-state  plow/soil  interaction  problem,  which  is  an  expensive  approach  from  a  computational 
point  of  view.  Further,  special  handling  techniques  are  necessary  to  the  development  of  the 
model  when  employing  commercial  software.  In  this  case,  these  techniques  are  mainly  due  to  the 
imposed  usage  of  a  high  order  finite  element  formulation  to  supporting  effective  stress  analysis. 
Data  management  is  impacted  by  large,  generated  files  due  to  midside  node  points.  Moreover, 
when  used  adjacent  to  a  sliding  interface,  a  well-known  anomaly  involving  high-order  element 
nodal  reaction  forces,  requires  the  addition  of  constrained  displacement  fields  on  element  faces 
which  are  in  contact  with  the  plow.  The  resulting  model  approach  is  feasible  but  technically 
constrained.  In  the  event  of  numerical  difficulty,  which  is  often  encountered  in  nonlinear  analysis 
computation,  there  are  few  options  to  the  proposed  approach  for  plow/soil  interaction  modeling 
when  using  commercial  finite  element  technology. 

A  summary  of  two  alternative  modeling  approaches  for  plow/soil  interaction  also 
developed  in  the  course  of  this  project  was  provided.  The  first  approach  was  a  semi-empirical 
model  calibrated  with  NFESC  laboratory  test  data,  and  it  produced  three-dimensional  results  in 
very  good  agreement  with  the  experimental  data  base.  The  other  approach  was  a  new  and  simpler 
rational  numerical  model  approach  that  promises  to  avoid  the  substantial  computational  cost  of 
multistep  transient  analysis  for  simulating  steady-state  plow/soil  interaction.  It  was  successfully 
demonstrated  using  a  two-dimensional  nonlinear  plane  strain  plow/soil  interaction  model  with  a 
simplified  viscoelastic  material  model  for  soil.  Steady-state  behavior  of  plow/soil  interaction  was 
captured  in  a  single  computational  finite  element  solution  step.  This  alternative  to  multistep 
models  requires  further  development  including  implementation  of  an  elasto-plastic  effective 
stress  soil  model  and  a  three-dimensional  finite  element  analysis  framework.  This  is  feasible 
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using  special  purpose  software  developed  jointly  at  the  University  of  California  at  Davis  and  at 
NFESC  under  ONR  sponsorship.  The  software  features  dynamic  effective  stress  analysis  with 
modem  soil  constitutive  laws,  low-order  finite  element  formulations  for  soil,  and  a  modular 
architecture  for  producing  and  testing  alternative  algorithms  and  achieving  numerical  solutions 
for  a  broad  range  of  structural/geotechnical  problems. 

The  technical  approach  to  plow/soil  interaction  modeling  (the  primary  focus  of  this 
study),  considered  a  two-dimensional  multistep  plow/soil  interaction  model  demonstration.  This 
model  featured  a  geometrically  nonlinear  plane  strain  formulation,  and  a  Drucker-Prager 
plasticity  model  for  soil  behavior.  Preliminary  studies,  which  employed  simple  constitutive 
models  for  soil,  were  initially  encouraging,  demonstrating  that  a  rigid  plow  blade  penetrating  into 
soil  at  constant  velocity  could  be  simulated  successfully.  Simultaneously  in  a  parallel 
development,  preliminary  calibration  of  a  Drucker-Prager  constitutive  soil  model  demonstrated 
that  experimental  triaxial  test  data  for  seafloor  soil  relating  effective  stress  to  soil  deformation 
could  be  reasonably  replicated. 

However,  replacing  the  simplified  soil  material  models  with  a  calibrated  Drucker-Prager 
effective  stress  soil  model  to  achieve  a  satisfactory  plow/soil  interaction  model  was  not 
straightforward.  Several  technical  difficulties  remained  and  had  to  be  surmounted  to  achieve  a 
reasonable  working  model  of  plow/soil  interaction.  These  technical  difficulties  included 
appropriate  modeling  of  soil  tearing  or  failure  at  the  tip  of  the  plow  blade,  allowance  for  overly 
constrained  drainage  path  lengths  for  the  porous  medium  in  a  two-dimensional  model,  and 
unavoidable  numerical  experimentation  with  solution  control  parameters  to  achieve  stable 
solutions  and  acceptable  compute  times. 

The  behavior  of  the  two-dimensional  plow/soil  interaction  model  was  initially  analyzed 
for  a  nominal  plow  velocity,  a  standard  blade  configuration,  and  a  frictionless  plow/soil  interface 
condition.  Selected  results  from  the  analysis  pertaining  to  both  the  mechanics  of  plow/soil 
interaction  and  the  veracity  of  the  model  included  the  following  observations: 

1 .  Displacement  of  soil  surrounding  the  plow  under  steady-state  conditions  revealed  that 
only  the  tip  of  the  plow  blade  is  in  contact  with  the  soil.  Further,  a  permanent  groove 
or  wake  is  formed  behind  the  plow  tip  which  is  in  qualitative  agreement  with 
experimental  observation  of  the  groove  produced  by  the  physical  model  in  the 
laboratory  apparatus.  The  model  strictly  predicts  that  this  space  is  void  of  soil 
including  its  pore  water.  However,  pore  water  under  reduced  pressure  may  reasonably 
be  assumed  to  occupy  the  groove  or  wake. 

2.  The  steady-state  tow  force  acting  on  a  physical  model  of  a  thin  blade  plow  moving  in 
a  given  soil  was  reliably  measured  in  the  laboratory,  but  cannot  be  computed  by  the 
plow/soil  interaction  model  due  to  its  two-dimensional  nature.  However,  unit  plow 
forces  were  predicted  which  were  useful  to  the  understanding  of  plow/soil  interaction, 
and  which  complemented  information  available  from  the  laboratory  apparatus.  For 
example,  the  model  provided  information  on  the  lateral  or  side  force  acting  on  the 
plow  which  cannot  be  measured  by  the  laboratory  apparatus. 

3.  Predicted  strain  response  in  the  soil  revealed  how  the  trench  is  formed  from  residual 
or  plastic  strain,  and  showed  that  the  soil  has  undergone  permanent  volumetric 


expansion,  as  should  be  anticipated  for  a  dilatant  seafloor  soil  subjected  to  shear 
stress.  Predicted  effective  stress  response  in  the  soil  showed  a  concentration  of  stress 
acting  at  the  shoulder  of  the  plow  tip  which  reveals  why  steel  blades,  after  continuous 
usage,  have  an  observed  distinct  pattern  of  wear  on  the  shoulder  of  the  plow  tip. 
Moreover,  the  computed  stress  response  showed  that  the  plow  tip  must  overcome  soil 
which  has  been  made  stronger  by  the  shearing  action  of  the  tip  itself,  just  ahead  of  the 
plow,  and  provided  an  explanation  of  why  high  tow  forces  observed  in  earlier  field 
operations  were  not  predicted. 

4.  In  comparison  to  strain  and  effective  stress,  both  of  which  exhibited  high  gradients, 
the  pore  pressure  gradient  response  was  very  diffuse.  This  tended  to  make  it  more 
difficult  to  simulate  the  pore  pressure  response  in  plow/soil  interaction  due  to  the 
presence  of  finite,  artificial  boundaries  in  the  model.  However,  the  predicted  pore 
pressure  response  in  the  region  of  interest  surrounding  the  plow  tip,  was  nonetheless 
reasonable  and  qualitatively  consistent  with  measurements  made  with  the  laboratory 
apparatus. 

Also,  a  subsequent  parameter  study  with  the  two-dimensional  plow/soil  interaction  model 
demonstrated  the  potential  engineering  utility  of  a  future  rational  three-dimensional  model.  The 
effects  of  plow  velocity,  interface  friction,  and  shape  of  the  plow  tip  were  investigated.  The 
results  are  presented  in  the  following: 

1.  Unit  tow  force  predicted  by  the  model  was  a  linear  function  of  plow  velocity,  in 
agreement  with  the  linear  relationship  between  tow  force  and  velocity  measured  with 
the  laboratory  apparatus.  However,  any  further  direct  comparison  of  the  data  is 
misleading  because  the  unit  tow  force  predicted  by  the  two-dimensional  model  and 
(total)  tow  force  measured  in  the  laboratory  are  different  quantities.  Further,  critical 
velocity,  as  observed  experimentally,  was  noted  in  the  range  of  velocities  studied  with 
the  model.  This  was  most  likely  due  to  two  reasons:  (a)  the  model  contained  no 
provision  for  cavitation  of  the  pore  fluid  and  (b)  a  proper  failure  condition  for  soil  at 
the  plow  tip  has  as  yet  to  be  incorporated  into  the  constitutive  model.  In  the 
neighborhood  of  the  plow  tip,  resisting  principal  stresses  in  the  soil  did  increase 
toward  an  apparent  bound,  while  maximum  shear  stress  decreased,  as  the  steady-state 
plow  velocity  increased.  The  decrease  in  pore  pressure  near  the  plow  tip  was  linear 
with  plow  velocity.  The  decrease  was  without  an  apparent  bound  in  the  absence  of  a 
cavitation  pressure  cutoff  in  the  model. 

2.  Unit  tow  force  varied  linearly  with  the  coefficient  of  friction  between  the  plow  blade 
and  soil.  For  example,  10  percent  friction  increased  the  tow  force  by  nearly  35  percent 
over  the  frictionless  case,  while  the  lateral  force  increased  only  very  modestly. 
Further,  10  percent  friction  caused  the  maximum  shear  stress  in  the  soil  near  the  plow 
tip  to  increase  dramatically  by  a  factor  of  3. 

3.  The  shape  of  pointed,  plow  blade  tip  profiles  is  clearly  an  important  engineering 
design  parameter.  However,  it  is  also  an  important  operational  parameter  since. 
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during  usage  in  the  field,  wear  may  be  anticipated  which  will  cause  its  shape  to 
change  with  corresponding  effects  on  performance  and  efficiency.  In  the  absence  of 
friction,  the  shape  parameter  study  found  that  the  unit  plow  tow  force  changed 
dramatically  with  blade  shape,  while  the  unit  plow  lateral  force  was  nearly 
independent  of  blade  shape.  The  unit  tow  force  increased  somewhat  faster  than 
linearly  with  the  size  of  the  included  angle  of  the  tip.  For  example,  it  was  halved  by 
halving  the  standard  60-degree  angle  of  the  plow  tip,  while  it  increased  75  percent  for 
a  50  percent  increase,  corresponding  to  a  dull  blade  tip  (90-degree  angle).  Further, 
because  the  maximum  lateral  contact  force  was  much  reduced  for  the  sharp  blade,  a 
sharper  blade  entry  would  potentially  sustain  much  less  shoulder  wear  on  the  tip. 
However,  the  dull  blade  produced  larger  shear  stress  and  larger  negative  pore  pressure 
in  the  soil.  At  the  same  towing  velocity,  it  would  be  potentially  more  effective  in  soil 
trenching  albeit  at  the  expense  of  higher  tow  forces. 

Ultimately,  a  three-dimensional  numerical  plow/soil  interaction  model  would  be 
necessary  to  accurately  predict  the  response  of  a  plow  moving  with  steady-state  velocity  through 
the  seafloor.  In  spite  of  the  high  technical  risk  in  the  development  of  such  a  model,  the  two- 
dimensional  model  presented  in  this  study  adequately  demonstrates  that  it  would  complement  the 
existing  knowledge  base  supporting  cable  burial  operations.  Often,  useful  scientific  and 
engineering  information  cannot  be  as  easily  obtained  by  alternative  methods  especially  when  the 
behavior  and  response  of  the  plow  depends  critically  on  either  site-specific  information  or  on 
particular  design  details  of  plow  blades. 
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7.0  RECOMMENDATIONS 


Future  cable  burial  technology  requirements  for  numerical  simulation  of  plow/soil 
interaction  should  be  addressed  by  focusing  on  the  development  of  three-dimensional  finite 
element  solutions.  If  both  transient  and  steady-state  solutions  are  required,  then  continued 
application  of  commercial  finite  element  technology  and  multistep  solution  approaches  together 
with  super  computer  resources  would  be  needed.  On  the  other  hand,  if  only  steady-state  solutions 
are  required,  then  the  recommended  approach  is  to  complete  the  development  of  the  proposed 
single-step  finite  element  solution  approach,  which  is  more  computationally  efficient. 

The  single-step  finite  element  approach  is  more  likely  to  result  in  an  effective  engineering 
analysis  and  design  tool  because  it  would  require  only  workstation  or  personal  computer 
resources.  Development  of  supporting  special-purpose  software  is  also  recommended  to  enhance 
that  goal.  Moreover,  this  approach  would  provide  for  application  of  classical  critical  state  soil 
models  and  modem  advanced  soil  constitutive  models,  such  as  the  bounding  surface  plasticity 
model  which  was  developed  by  the  Navy  and  which  is  generally  unavailable  with  commercial 
software. 

It  is  also  recommended  that  additional  experiments  be  conducted  to  determine  the  failure 
mechanism  of  soil  at  the  tip  of  the  plow.  This  mechanism  has  not  been  given  sufficient  attention, 
though  it  is  important  to  the  material  modeling  process  irrespective  of  the  numerical  modeling 
approach.  This  is  especially  important  to  establishment  of  the  upper  bound  in  plow  resistance 
which  governs  actual  plow  design.  Also,  further  experiments  specifically  designed  for  validating 
models  of  plow/soil  interaction  are  recommended. 

Completion  of  the  single-step  finite  element  solution  approach  for  steady-state  plow/soil 
interaction  problems  requires  the  following  recommended  steps  be  taken: 

1.  Perform  additional  experiments  with  the  laboratory  apparatus,  along  with  further 
instrumentation,  to  define  failure  modes  of  soil  at  the  tip  of  rigid  penetrators,  and  to 
provide  additional  validation  data  for  plow/soil  interaction  models. 

2.  Develop  failure  models  from  these  experiments  and  implement  them  into  the  single- 
step  finite  element  model. 

3.  Replace  the  existing  simplified  viscoelastic  soil  model  with  the  bounding  surface 
plasticity  soil  model  in  the  single-step  finite  element  software.  A  critical  state  soil 
model  has  in  fact  already  been  implemented  for  undrained  conditions  in  this  manner. 

4.  Extend  the  solution  framework  of  the  single-step  finite  element  software  from  two 
dimensions  to  three  dimensions. 

5.  Perform  related  numerical  studies  to  verify  implementation  of  algorithms,  and 
validate  the  single-step  finite  element  model  using  experimental  data  from  the 
laboratory  apparatus. 
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A  related  requirement  for  simulating  two-dimensional  axisymmetric  penetration  response 
of  penetrometers  in  saturated  sand  is  similar  to  the  requirement  addressed  in  this  study  for 
simulating  plow/soil  interaction.  Moreover,  a  two-dimensional  axisymmetric  formulation  of  the 
penetrometer  problem  is  both  natural  and  reasonable.  A  typical  impact  velocity  of  a  penetrometer 
is  50  fps,  sufficiently  low  that  the  problem  is  amenable  to  a  transient  solution  with  an  implicit 
finite  element  program.  This  has  been  the  approach  used  in  the  present  study  of  steady-state 
plow/soil  interaction  at  velocities  from  1  to  13  fps.  An  explicit  finite  element  solution  could  also 
be  recommended  for  the  transient  penetrometer  problem,  but  cannot  also  be  easily  recommended 
for  the  steady-state  problem. 

Therefore,  both  problems  are  amenable  to  dynamic  solutions  by  implicit  finite  element 
programs,  both  involve  seafloor  soil,  and  both  require  modeling  of  frictional  interfaces  using 
slide  lines.  Consequently,  much  of  the  expertise  gained  in  the  study  of  plow/soil  interaction 
modeling  is  immediately  applicable  to  numerical  modeling  of  penetrometers.  A  general 
numerical  modeling  methodology  can  be  designed  to  address  both  problems  using  two-  and 
three-dimensional,  transient  and  steady-state,  nonlinear  dynamic  finite  element  technology. 

A  computer-based  design  and  analysis  tool  supporting  future  requirements  in  cable  burial 
technology  and  related  seafloor  geotechnical  operations  is  recommended.  This  recommendation 
will  have  to  be  considered  and  proposed  in  more  detail,  but  sufficient  engineering  resources  (i.e., 
experimental,  theoretical,  and  numerical  modeling  expertise)  are  in  place  to  realize  this  product 
within  a  reasonable  time  frame  and  for  a  reasonable  research  and  development  investment. 
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Appendix  A 

DISCUSSION  OF  SOME  TECHNICAL  ISSUES  AND  PROJECT  CONSTRAINTS 

The  impact  of  limited  project  funding  and  a  compressed  project  schedule  for  determining 
the  capability  of  commercial  finite  element  technology,  as  represented  by  the  ABAQUS 
program,  to  model  plow/soil  interaction  is  discussed.  The  discussion  is  organized  into  six 
sections:  (1)  Material  Type  and  Soil  Behavior,  (2)  Material  Model,  (3)  Finite  Element  Model, 
(4)  Parameter  Study,  (5)  U.  C.  Davis  Studies,  and  (6)  Contingency  Model  Study.  A  summary  of 
these  issues  and  their  effect  on  evaluation  of  numerical  modeling  with  ABAQUS  is  provided. 

1.1  Material  Type  and  Soil  Behavior 

The  plow  is  required  to  traverse  medium  and  fine  seafloor  sands  under  the  following 
conditions:  100  percent  saturation,  high  density,  drained,  partially  drained,  and  undrained 
conditions.  Emphasis  will  be  given  to  100  percent  saturated  medium  sand  under  partially  drained 
conditions. 

Required  soil  behavior  necessitates  a  two-phase,  effective  stress  soil  model  approach 
prior  to  cavitation,  and  possibly  a  three-phase  system  approach  if  cavitation  occurs.  The  material 
is  subjected  to  both  deviatoric  and  dilatation  stresses,  as  well  as  pore  pressure  changes. 

Cavitation  cannot  be  directly  addressed  in  this  study  because  ABAQUS  does  not  provide 
for  it.  Thus,  the  analysis  remains  limited  to  a  two-phase,  highly  dense  soil,  with  100  percent 
saturation  conditions.  Flow  conditions  are  discussed  below  in  connection  with  the  material 
model.  Sands  are  frictional  materials.  Therefore,  a  nonassociative  plasticity  model  will  be 
addressed.  Associative  models  will  not  be  studied. 

1.2  Material  Model 

An  assessment  of  commercial  capability  could  be  incomplete  if  some  applicable 
capabilities  are  excluded  from  the  study.  In  the  case  of  ABAQUS,  three  soil  models  for  sands  are 
available:  (1)  Modified  Drucker-Prager  Model,  (2)  Modified  Cap  Model  and  (3)  Modified  Cam- 
Clay  Model. 

The  Modified  Cap  Model  includes  three  failure  surfaces,  a  compaction/dilation  limitation 
surface  to  model  volume  change,  a  curved  Drucker-Prager  surface  so  that  the  direction  of 
material  loading  is  permitted  to  vary,  and  a  transition  surface  between  the  latter  two  surfaces 
which  allows  for  modeling  material  softening.  Most  of  the  calibration  effort  is  related  to  the 
Drucker-Prager  shear  surface.  Little  more  effort  is  needed  beyond  that  required  to  calibrate  the 
Modified  Drucker-Prager  Model.  Though  study  of  the  Modified  Cap  Model  would  be  useful,  it  is 
also  true  that  available  material  calibration  data  is  limited.  The  Modified  Cam-Clay  Model 
includes  an  elliptical  failure  surface  (in  the  meridian  plane).  It  also  includes  a  dependence  on  the 
third  stress  invariant.  However,  this  model  is  believed  to  do  essentially  the  same  things  as  the  cap 
model  and  its  inclusion  in  the  study  would  be  somewhat  duplicative. 

Only  one  applicable  material  model,  the  Modified  Drucker-Prager  will  be  addressed.  The 
numerical  model  will  be  limited  to  the  extended  Drucker-Prager  plasticity  model  for  the  soil 
skeleton,  along  with  the  pore  water  pressure  variable  formulation.  In  the  drained  condition,  the 
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model  is  employed  solely,  without  pore  water  pressure  variables.  In  the  other  two  conditions,  the 
material  model  is  employed  along  with  the  pore  water  pressure  variables.  The  partially  drained 
condition  can  be  handled  naturally  through  specification  of  the  boundary  value  problem  in  which 
pore  water  pressure  is  modified  by  whatever  dissipation  mechanisms  are  present.  The  undrained 
condition  is  specified  by  an  arbitrarily  low  permeability  coefficient.  In  this  regard  there  is  limited 
impact  due  to  time  constraints  in  the  event  of  technical  problems  involving  some  of  these 
conditions.  For  example,  the  undrained  condition  may  be  dropped  from  the  task  if  numerical 
problems  involving  oscillation  in  pore  water  pressures  arise. 

Only  a  nonassociative  flow  rule  will  be  considered  because  of  the  material’s  frictional 
nature.  The  computations  are  more  lengthy  with  a  nonassociative  flow  rule  because  the  resulting 
equations  of  equilibrium  are  nonsymmetric.  An  associative  flow  rule  would  in  effect  discount 
volume  change  in  the  soil  which  is  otherwise  believed  to  be  important  in  the  present  problem. 

1.3  Finite  Element  Model 

It  is  believed  that  a  three-dimensional  stress  state  prevails  in  the  vicinity  of  the  plow 
which  cannot  be  captured  by  a  two-dimensional  model  without  potentially  far  reaching 
assumptions  on  the  applicability  of  the  results.  It  is  also  not  clear  that  friction  on  the  plow  may  be 
safely  neglected.  Friction  should  be  retained  in  the  study.  Yet,  it  is  known  that  this  complicates 
the  modeling  process.  Further,  large  deformations  clearly  exist  in  the  neighborhood  of  the  plow. 
Also,  inertia  forces  due  to  soil  particle  accelerations  are  potentially  an  issue,  for  if  they  are 
significant  they  can  affect  pore  pressure,  effective  stress  and  resistance  force  on  the  plow. 

A  two-dimensional  plane  strain  plow/soil  interaction  model  is  not  easily  defended. 
ABAQUS  does  not  support  plane  stress  formulations  for  soil  consolidation  problems.  A  three- 
dimensional  model  simply  cannot  be  considered  beyond  some  initial  concepts  due  to  time  and 
funding  constraints.  Therefore,  emphasis  is  given  to  a  two-dimensional  plane  strain  finite 
element  model  of  plow/soil  interaction. 

An  attempt  will  be  made  to  retain  friction  in  the  study  though  numerical  difficulty  is  a 
distinct  possibility  when  it  is  included.  In  ABAQUS  there  are  finite  element  families  for 
addressing  node-to-node  and  surface-to-surface  interaction  under  small  or  finite  sliding 
assumptions  with  or  without  friction.  Emphasis  will  be  given  to  a  frictionless  plow/soil 
interaction  model. 

Large  deformations  will  be  included  in  conjunction  with  the  Drucker-Prager  plasticity 
model,  so  no  scope  changes  are  anticipated  in  this  respect.  However,  an  inconsistency  may  exist 
due  to  specification  of  large  deformations  in  conjunction  with  a  material  model  based  upon  small 
strain  measures.  It  is  believed  that  this  approach  is  nonetheless  worthwhile,  for  the  alternative  is 
to  limit  the  analysis  to  small  deformations  which  is  incorrect  to  first  order. 

Quasi-static  transient  consolidation  analysis  will  be  pursued.  Dynamic  consolidation 
analysis  cannot  be  evaluated  since  inertia  forces  are  neglected  in  ABAQUS  for  coupled  effective 
stress/porous  media  flow  problems. 

1.4  Parameter  Study 

There  are  many  parameters  to  be  considered  in  seeking  to  establish  adequacy  of  the 
model  and  to  learn  from  the  model.  Limited  resources  coupled  with  very  large  run  times  and 
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heavy  post-processing  and  data  analysis  burdens  will  allow  only  so  many  parameters  to 
evaluated.  An  attempt  will  be  made  to  include  limited  parameter  studies  on  the  effects  of  plow 
resistance  from  steady-state  velocity,  friction,  blade  sharpness,  and  permeability. 

1.5  U.  C.  Davis  Studies . 

Two  studies  were  commissioned  early  in  the  numerical  model  project  at  the  University  of 
California  at  Davis.  The  main  task  of  these  studies  was  to  provide  discussions  of  the  relevant  soil 
mechanics  for  plow/soil  interaction  and  numerical  modeling  issues  to  guide  the  present 
development  of  a  numerical  model.  If  the  numerical  model  project  schedule  were  to  have  been 
accelerated,  the  main  purpose  of  these  studies  would  have  been  thwarted.  However,  the  project 
could  not  be  accelerated  as  results  were  not  forthcoming  early  enough.  Nonetheless,  these  two 
studies  were  eventually  published  (Kutter  and  Voss,  1995  and  Herrmann  and  Mello,  1994)  and 
provided  useful  guidance  as  initially  planned’. 

1.6  Contingency  Model  Study 

Given  the  importance  of  the  problem  to  the  Navy  and  the  early  uncertainties  in  modeling 
success,  contingency  modeling  plans  were  to  be  made  and  conducted  somewhat  in  parallel  with 
primary  numerical  modeling  effort.  However,  project  constraints  largely  eliminated  real  progress 
on  these  contingencies.  Further,  success  of  the  model  remains  a  subjective  issue  in  the  absence  of 
specific  criteria.  Generally,  it  has  meant  the  ability  of  the  model  to  replicate  measured  response 
from  the  laboratory  apparatus. 

If  the  two-dimensional  plow/soil  interaction  model  fails  to  agree  with  experimental 
observation  and  measurement  and  a  three-dimensional  model  is  deemed  necessary  as  a  result, 
then  the  course  of  action  would  be  to  complete  the  originally  proposed  study  which  called  for 
construction  of  a  three-dimensional  model  and  substantial  usage  of  super  computer  resources.  If 
the  two-dimensional  model  fails  for  other  reasons,  then  the  option  would  be  to  consider  possibly 
two  other  courses  of  action  depending  on  the  difficulty.  First,  if  the  material  models  in  ABAQUS 
do  not  suffice,  then  a  modem  material  model  can  be  implemented  in  ABAQUS  so  that  other 
important  modeling  resources  of  ABAQUS  may  still  be  exploited.  If  the  contact  algorithms  and 
the  provisions  for  large  sliding  interfaces,  etc.,  contained  in  ABAQUS  proved  inadequate,  then 
an  approach  using  Geneseos,  an  in-house  nonlinear  finite  element  program,  can  be  considered  in 
place  of  commercial  technology. 

Three  noncommercial  finite  element  technology  alternatives  to  the  main  numerical  model 
approach  (featuring  commercial  technology),  were  considered  as  contingency  model  studies: 

1.  Application  of  ABAQUS  with  the  U.  C.  Davis  bounding  surface  plasticity  soil  model 
replacing  the  Drucker-Prager  plasticity  soil  model.  Quasi-static  transient 
consolidation  analysis  could  then  be  resumed. 

2.  Application  and  further  development  of  Geneseos,  which  has  the  bounding  surface 
plasticity  soil  model,  but  which  requires  large  deformation  and  suitable  interface 
modeling  capabilities.  Fully  dynamic  transient  soil  consolidation  analysis  could  then 
be  pursued. 
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3.  Application  of  U.C.  Davis’  steady-state  finite  element  analysis  formulation  of  the 
plow  problem  (Herrmann  and  Mello,  1994).  Fully  dynamic  steady-state  consolidation 
analysis  could  then  be  pursued. 

The  current  numerical  model  study  incorporates  the  Drucker-Prager  model  to  represent 
the  response  of  the  soil.  If  this  model  is  deemed  inadequate,  the  CAP  plasticity  model  must  also 
be  considered  to  fully  evaluate  the  current  approach  of  using  commercial  technology.  If  both  of 
these  soil  models  prove  inadequate  and  if  the  U.C.  Davis  bounding  surface  plasticity  can 
potentially  overcome  their  deficiencies,  then  it  could  be  implemented  with  ABAQUS  via  a 
UMAT  procedure  (a  user-written  subroutine  for  incorporating  user  supplied  material  laws).  Such 
a  procedure  has  in  the  past  been  written  for  an  earlier  version  of  the  bounding  surface  plasticity 
soil  model. 

The  current  quasi-static  transient  consolidation  analysis  approach  neglects  inertia  effects. 
If  these  effects  are  determined  to  be  significant,  ABAQUS  would  have  to  be  abandoned  for  it  has 
no  dynamic  soil  consolidation  capability.  An  alternative  course  of  action  would  be  additional 
development  of  Geneseos  to  support  the  present  problem.  It  has  the  bounding  surface  plasticity 
soil  model,  corresponding  low  order  two-dimensional  soil  elements  and  either  quasi-static  or 
fully  dynamic  consolidation  analysis  capability.  The  three  main  development  tasks  for  this 
approach  would:  (1)  add  a  slide-line  interface  capability,  (2)  add  large  deformation  capability, 
and  (3)  update  existing  three-dimensional  elements  to  soil  elements  by  providing  for  pore 
pressures. 

All  of  the  numerical  model  approaches  considered  to  this  point  occur  in  the  time  domain. 
Accordingly,  the  plow  must  enter  the  soil  and  be  towed  through  a  transient  phase  until  steady- 
state  forces  are  established.  An  alternative  steady-state  finite  element  approach  is  under 
development  and  could  be  pursued  if:  (1)  the  time  domain  approaches  fail  to  converge  on  a 
solution,  or  (2)  the  time  domain  approaches  converge  but  are  computationally  prohibitive  and 
prevent  parameter  studies  in  support  of  engineering  design  to  be  conducted. 

This  approach  also  requires  substantial  software  development,  but  it  also  provides  the 
greatest  potential  return.  If  successful,  a  nonlinear  steady-state  analysis  of  the  plow  traversing  the 
soil  could  be  as  inexpensive  as  one  time  step  in  a  nonlinear  dynamic  transient  analysis.  Thus 
parameter  studies  to  examine  the  effects  of  numerical  and  physical  variables  would  be  supported. 
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Appendix  B 


STATEMENT  OF  WORK  FOR  U.  C.  DAVIS  STUDIES 


1.0  INTRODUCTION  AND  BACKGROUND 

1.1  The  Naval  Civil  Engineering  Laboratory  (NCEL)  within  L40  has  been  tasked  to 
investigate  high  performance,  lightweight  burial  plows  suitable  for  burial  of  small  diameter 
cables  in  specific  seafloor  soils.  In  conjunction  with  this,  L51  has  been  tasked  to  develop  a 
numerical  model  of  the  plow/soil  interaction  process,  the  goal  of  which  is  to  calculate  the  soil 
resisting  force  on  the  blade  of  the  plow  during  towing  at  a  constant  velocity  (0  to  10  knots 
range). 

1.2  The  plow/soil  interaction  process  is  physically  complicated  and  requires  geotechnical 
clarification.  A  soil  mechanics  study  is  required  which  will  include  a  discussion  and  description 
of  the  relevant  physics  involved  including  the  identification  of  parameters  and  their  potential 
effects,  a  discussion  of  the  relevance  of  existing  bounding  surface  soil  models  for  clay  and  sand 
materials,  a  discussion  and  description  of  the  requirements  for  a  valid  numerical  model  of 
plow/soil  interaction,  and  a  discussion  of  the  potential  for  using  existing  computer  programs  for 
the  numerical  model. 


2.0  SCOPE 

2.1  The  contractor  shall  provide  all  labor  and  materials  required  to  prepare  a  discussion  on 
the  relevant  physics  of  plow/soil  interaction,  necessary  modifications  to  NCEL/UCD  bounding 
surface  soil  models,  requirements  for  numerical  modeling,  and  potential  for  using  commercial  or 
available  computational  mechanics  software  (ABAQUS,  DYNA3D,  Geneseos,  DYSAC,  etc.). 

2.2  The  key  element  of  the  numerical  model  is  the  constitutive  model  used  for  describing  the 
stress-strain  behavior  of  the  seafloor  soil.  Since  NCEL  is  tasked  to  develop  this  model,  the 
proposed  study  must  be  conducted  within  the  framework  of  NCEL’ s  soil  models.  These  models 
have  been  developed  in  conjunction  with  Professor  Kutter  and  his  colleagues  at  UCD.  (Refer  to 
the  recently  completed  contract  N47408-89-C-1058,  Validation  of  a  Proposed  Rational  Material 
Model  Characterization  of  Granular  Soils.)  Their  data  requirements,  performance  capabilities, 
and  required  modifications  shall  be  considered  while  forming  recommendations  for  the 
numerical  modeling  task  at  NCEL. 


3.0  GOVERNMENT  FURNISHED  INFORMATION 

3.1  Plow/Soil  Interaction  Data.  The  government  will  provide  plow/soil  interaction 
experimental  data  and  relevant  technical  literature  as  guidance  to  the  contractor  for  preparing  his 
study. 
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4.0  TECHNICAL  TASKS 


4.1  The  study  shall  include  a  discussion  and  description  of  the  relevant  physics  and  soil 
mechanics  behavior  involved  in  the  plow/soil  interaction  process.  The  study  shall  include  the 
identification  of  key  parameters  and  a  discussion  of  their  importance.  Parameters  and  variables 
such  as  soil  type,  particle  size,  permeability,  relative  density,  pore  pressure,  cavitation  and  shear 
strength  shall  be  included.  This  task  shall  form  the  main  part  of  the  study. 

4.2  The  study  shall  be  conducted  within  the  framework  of  the  NCEL/UCD  bounding  surface 
plasticity  soil  model  for  clay  and  the  bounding  surface  hypoplasticity  soil  model  for  sand.  The 
data  requirements,  performance  capabilities  and  any  required  modifications  for  these  models 
shall  be  considered  for  support  of  a  plow/soil  numerical  modeling  task  at  NCEL.  Soil 
constitutive  model  recommendations  shall  be  given. 

4.3  The  study  shall  include  a  discussion  and  description  of  the  requirements  for  a  valid  finite 
element  numerical  model  of  plow/soil  interaction.  In  addition  to  the  requirements  of  the 
constitutive  model,  other  finite  element  modeling  considerations  shall  be  considered  such  as 
transient  loading,  steady-state  loading,  boundary  conditions,  slide  lines,  small  and  large 
deformations,  etc.  Approaches  other  than  finite  element  modeling  may  also  be  included  in  a 
review  of  alternative  approaches  to  plow/soil  interaction.  Finite  element  modeling 
recommendations  shall  be  provided. 

4.4  The  study  shall  include  an  evaluation  of  the  potential  for  using  commercial  and  available 
general  purpose  finite  element  computer  programs  such  as  ABAQUS  and  DYNA3D,  and 
available  research/developing  programs  such  as  Geneseos  (NCEL)  and  DYSAC  (UCD)  for 
modeling  the  plow/soil  interaction  problem.  Computer  program  recommendations  shall  be  given. 


5.0  DELIVERABLES 

5.1  Technical  Report  on  Plow/Soil  Interaction  and  Related  Numerical  Modeling  Issues. 

A  final  technical  report  on  the  relevant  soil  mechanics  issues  involved  with  plow/soil  interaction 
shall  be  prepared.  The  study  shall  also  include  recommendations  on  constitutive  models, 
numerical  modeling,  and  finite  element  computer  programs.  The  report  shall  include  an  abstract 
of  the  findings. 

5.2  Six  (6)  copies  of  the  camera-ready  final  technical  report  shall  be  provided. 

6.0  SCHEDULE 

6.1  The  final  technical  report  shall  be  delivered  within  eight  (8)  months  of  the  start  of  the 
contract. 
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Appendix  C 


ANNOTATED  INPUT  DATA  FILE  FOR  PLOW/SOIL  INTERACTION  MODEL 


A  description  of  a  typical  input  data  file  is  provided  for  documentation  of  the  plow/soil 
interaction  model  and  in  the  event  that  future  simulations  are  required.  In  the  following  input  file 
listing,  any  command  preceded  by  **  is  ignored  by  ABAQUS,  i.e.,  “commented  out.”  The 
annotations  provided  are  italicized  and  enclosed  in  parentheses  following  each  command. 
Otherwise  the  command  lines  or  “cards”  are  exactly  as  they  appear  in  a  typical  input  data  file. 

^HEADING,  UNSYMM  (specifies  an  unsymmetrical  equation  solver  algorithm 

due  to  flow  rule) 

PLOW/SOIL  INTERACTION:  1 /2-inch  (this  heading  is  printed  on  output) 

pointed  blade  at  5  fps,  2-D  plane  strain 

model,  CPE8P  elements,  interface  friction 

coeff.=0.1,  Drucker-Prager  soil  w/o 

hardening/softening,  yield=14  psi 

compression,  with  nonassociative  flow 

coupled  pore  pressure/effective  stress 

analysis,  reduced  fluid  bulk  modulus  to 

compensate  for  no  cavitation  model. 


*NODE 

1,0.0,0.0 

145,0.0,9.0 

151,0.0,12.0 

3201,1.0,0.0 

3345,1.0,9.0 

3351,1.0,12.0 

4801,2.0,0.0 

4945,2.0,9.0 

4951,2.0,12.0 

5601,4.0,0.0 

5745,4.0,9.0 

5751,4.0,12.0 


(key  node  points  and  their  coordinates  in  inches  that 
define  a  rectangular  mesh) 

(node  1  is  the  origin  of  x-y-z  global  axes  which  is  placed 
at  the  start  end  of  mesh) 

(145  nodes,  including  midpoint  nodes,  to  span  first  9 
inches  along  positive  y-axis) 

(coarser  layout  of  6  nodes  from  9  -12  inches,  to  stop  end 
of  the  mesh) 

(specifies  a  duplication  of  the  above  node  pattern  across 
to  the  x=l  line) 

(repeat  of  the  second  card  of  the  list  for  nodes  along  the 
x=l  line  of  nodes) 

(repeat  of  the  third  card,  along  the  x=l  line  of  nodes) 
(specifies  a  coarser  mesh  between  x=l  and  x=2  inches) 
(repeat  of  the  second  card for  the  x-2  line  of  nodes) 
(repeat  of  the  third  card,  for  x=2  line  of  nodes) 

(specifies  an  even  coarser  mesh  between  x=2  and  x=4 
inches) 

(repeat  of  the  second  card  for  the  x=4  line  of  nodes) 
(repeat  of  the  third  card  for  the  x=4  line  of  nodes) 
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6001,6.0,0.0 

6145,6.0,9.0 

6151,6.0,12.0 

*NODE,NSET=PLOWNODE 

10001,0.0,-3.25 

*NODE,NSET=SYMMREST 

10002,0.0,6.0 

*NGEN 

1.145.1 
*NGEN 

145.151.1 

*NSET, GENERATE, NSET=CL 

1.151.1 
*NGEN 

3201.3345.1 
*NGEN 

3345.3351.1 

*NSET, GENERATE, NSET=T1 

3201.3351.1 
*NGEN 

4801.4945.1 
*NGEN 

4945.4951.1 

*NSET, GENERATE, NSET=T2 

4801.4951.1 
*NGEN 

5601.5745.1 
*NGEN 

5745.5751.1 

*NSET, GENERATE, NSET=T3 

5601.5751.1 
*NGEN 

6001.6145.1 
*NGEN 

6145.6151.1 

*NSET,GENERATE,N  SET=EDGE 

6001.6151.1 


(the  lower  right  corner  of  the  mesh  at  x=6  inches) 

(repeat  of  the  second  card  for  the  x=6  line  of  nodes) 
(repeat  of  the  third  card  to  the  stop  end  of  the  mesh) 

(a  key  node;  the  reaction  force  at  this  node  is  the  tow 
force) 

(defines  the  plow  node  at  the  centroid  of  the  rigid  body 
plow) 

(gives  name  to  a  node  on  the  rigid  surface  along  the 
centerline) 

(defines  a  node  at  the  center  of  the  12-inch  rigid  surface 
along  the  centerline) 

(generates  the  nodes  between  nodes  1  and  145, 
incrementing  by  1) 

(generates  the  nodes  between  145  and  151,  incrementing 

byl) 

(names  the  set  of  nodes  1  - 151  on  the  centerline,  CL) 

(generates  the  nodes  between  3201  and  3345, 

incrementing  by  1) 

(generates  the  nodes  between  3345  and  3351, 

incrementing  by  1) 

(names  the  set  of  nodes  3201  -  3351,  Tl) 

(generates  the  nodes  between  4801  and  4945) 

(generates  the  nodes  between  4945  and  4951) 

(names  the  set  of  nodes  4801  and  4951,  T2) 

(generates  the  nodes  between  5601  and  5745) 

(generates  the  nodes  between  5745  and  5751) 

(names  the  set  of  nodes  5601  and  5751,  T3) 

(generates  the  nodes  between  6001  and  6145) 

(generates  the  nodes  between  61 45  and  61 51) 

(names  the  set  of  nodes  6001  and  6151,  EDGE) 
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*NFILL,N  SET=E  1 
CL, Tl, 16, 200 
*NFILL,NSET=E2 
T1,T2, 8,200 
*NFILL,NSET=E3 
T2,T3 ,4,200 
*NFILL,NSET=E4 
T3,  EDGE,  2,200 

*NSET,GENERATE,NSET=ALL 

1,6151,1 

*N  SET, GENERATE, NSET^START 
1,6001,200 

*N  SET,GENERATE,NSET=STOP 
151,6151,200 

*ELEMENT,TYPE=CPE8P 

1,1,401,403,3,201,402,203,2 
*  ELGEN,ELSET=S  OIL 


(names  the  set  of  nodes  between  sets  CL  and  Tl,  El) 

(names  the  set  of  nodes  between  sets  Tl  and  T2,  E2) 

(names  the  set  of  nodes  between  sets  T2  and  T2,  E3) 

(names  the  set  of  nodes  between  sets  T3  and  EDGE,  E4) 

(names  the  set  of  all  nodes,  ALL) 

(names  the  set  of  nodes  along  the  start  end,  START) 

(names  the  set  of  nodes  along  the  stop  end,  STOP) 

(specifies  an  eight-node  plane  strain  /pore  pressure 
element) 

(defines  the  connectivity  or  node  pattern  for  element  #1) 
(generates  elements  1  to  75  along  y-axis,  and  names  this 
set  SOIL 


1.75.2.1.15.400.100 
**ELSET,ELSET=ECL,  GENERATE 
**1,75,1 

*ELSET,ELSET=STOP,GEN 

75.1475.100 

*ELSET,ELSET=RIGHT,GEN 

1401,1475,1 

*  SOLID 

SECTION, ELSET=SOIL,MATERIAL=DP  1 . 

*  MATERIAL, NAME=  DPI 
*ELASTIC 

6000.0,0.3 

*DRUCKER  PRAGER 
46.2,1. ,21. 5 

*YIELD,TYPE=COMPRESSION 
**0.1, 0.0 
**14.0,0.005 
**21.0,0.01 
**25.0,0.015 


(names  centerline  elements  1  -  75,  set  ECL) 

(names  element  set  75  -1475  along  the  stop  end,  STOP) 

(names  set  of  elements  1401  -  1475,  RIGHT) 

(assigns  mat  7  property  set  DPI  to  elset  SOIL) 

(begin  the  definition  of  a  material  named  DPI) 

(specify  elastic  properties) 

(modulus  of  elasticity  =  6,000  psi;  Poisson  ratio  =  0.3) 
(specify  extended  Drucker  Prager  plasticity  model 
w/nonassociative  flow) 

(friction  angle  or  strength  =  46.2  deg.;  angle  of  dilation  = 
21.5  deg.) 

(specifies  deviatoric  yield  behavior  in  compression  for  DP 
mat  ’l) 

(table  of  values  for  deviatoric  stress  and  corresponding 
axial  plastic  strain, 

i.e.,  14  psi  @  0.005  plastic  strain,  etc.;  data  obtained 
from  triaxial  test  data; 

calibration  includes  strain  hardening  and  softening 
behavior; 

ABAQUS  could  not  handle  this  data  due  to  DP  singularity 
condition  in  tension; 
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**27.0,0.021 

**15.0,0.14 
14., 0.0 

*  PERMEABILITY 

0.011,0.5 

0.015,0.563 

0.019,0.6 


strain  hardening/softening  was  discarded  ,  only  residual 

stress  behavior  retained; 

residual  stress  apparent  in  triaxial  data) 

(elastic  perfectly  plastic  behavior  specified,  residual 
stress  = 14  psi  for  all  strain) 

(specifies  soil  permeability  coefficient  as  a  function  of 
void  ratio) 

(table  of  values  gives  coefficient  and  void  ratio) 

(for  example,  k  —  0.015  and e  =  0.563) 

(this  is  measured  data  idealized  by  following  relationship: 

k  =  0.087  e^ ,  with  units  of  in/s  for  k) 


0.030,0.7 


0.044,0.8 

0.063,0.9 

0.087,1.0 

*POROUS  BULK  MODULI 
0.0,305000.0 

*  INITIAL  CONDITIONS, TYPE=RATIO 
ALL, 0.563 


(specifies  a  finite  bulk  modulus  for  soil  phases) 

(0. 0  implies  infinite  value  for  solid  phase  ;  bulk  modulus 
=  305,000 psi  for  fluid) 

(specifies  initial  void  ratio) 

(node  points  in  NSET  ALL  are  assigned  an  initial  value  of 
0.563) 


*  INITIAL  CONDITIONS, TYPE=PORE  (specifies  initial  pore  pressure) 


PRESSURE 
ALL, 0.0 

*  ELEMENT, TYPE=IRS2 1 

10001,1,3,10001 

*ELGEN,ELSET=TUNNEL 

10001,75,2,1 

*INTERFACE,ELSET=TUNNEL 

1.0 

*FRICTION,TAUMAX=l  000 
0.1 

*RIGID 

SURF  ACE, ELSET=TUNNEL,TYPE=SEGM 
ENTS,SMOOTH=0. 1 25 
START, -0.0025, 0.0044337567 

LINE,0.25, -0.44337567 
LINE, 0.25, -6.56 


(node  points  in  NSET  ALL  are  assigned  an  initial  value  of 
0. 0  psi) 

(planar  rigid  surface  element  for  use  with  first  order 
elements) 

(ele  #  10001;  node  1;  node  3,  rigid  body  reference  node 

10001) 

(generates  a  set  of  75  ele ’s  beginning  with  ele  # 1 0001) 
(master  ele  #10001;  no.  of  ele’s;  node  increment;  element 
increment) 

(defines  ele ’s  in  set  TUNNEL  to  be  contact  elements) 
(dummy  mechanical  property  for  contact  element  set) 
(Coulomb  friction  model  with  upper  bound  shear  stress 
=1000psi) 

(coefficient  of friction  between  rigid  plow  and  soil) 

(the  ele  set  TUNNEL 

is  to  be  used  in  conjunction  with  blade  rigid  surface 
defined  with 

tip  at  origin,  60  deg.  bevel  and  0.25  inch  half  thickness) 
(blade  length  is  6.25  inches ) 
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LINE,0.0,-6.56 
*ELEMENT,TYPE=IRS2 1 

10101,3,5,10002 

*NSET,NSET-FIRST,GENERATE 

1.149.2 

*  NS  ET,NS  ET=L  AS  T, GENERATE 

3.151.2 

*NSET,NSET=MID, GENERATE 

2.150.2 

*MPC 

LINEAR, MID, FIRST,LAST 

*ELGEN,ELSET=SYMMREST 

10101,74,2,1 

*ELSET,ELSET=START,GEN 

1,1401,100 

*INTERFACE,ELSET=SYMMREST 

1.0 

**  SURFACE  CONTACT,SOFTENED 
**0.0625,5000 

*BOND  SURFACE 

*RIGID 

SURFACE, ELSET=SYMMREST,TYPE=SE 

GMENTS 

START,0.0,13.0 

LINE, 0.0, -1.0 

*NSET,  NSET=CRACKREF 
1 

* BOUNDARY 

STOP,2 

EDGE,1 


(rigid  surface  specification  closes  on  centerline) 

(planar  rigid  surface  element  for  use  with  first  order 
elements) 

(ele  #  10101;  node  3;  node  5,  rigid  body  reference  node 

10002) 

(set  of  bottom  corner  nodes  along  centerline  is  called 
FIRST) 

( node  1;  149  such  nodes;  node  increment  for  generating 
set  is  2) 

(set  of  top  corner  nodes  along  centerline  is  called  LAST) 
(node  3;  151  such  nodes;  node  increment  for  generating 
set  is  2) 

(set  of  midside  nodes  along  centerline  is  called  MID) 

(node  2;  150  such  nodes;  node  increment  for  generating 
set  is  2) 

(specifies  multi-point  constraint  for  sets  MID,  FIRST, 
LAST) 

(midside  nodes  interpolate  displ ’s  linearly  between  corner 
nodes) 

(generates  a  set  of  74  ele’s  beginning  with  ele  #  10101 
which  is  tied  to  ele  #  2) 

(master  ele  #  10101;  no.  of  ele’s;  node  increment;  element 
increment) 

(defines  14  elements  along  starting  edge  of  soil  region) 
(ele  # 1 ;  ele#1401,  ele  increment  =  100) 

(defines  ele ’s  in  SYMMREST  to  be  contact  elements) 
(dummy  mechanical  property  for  contact  element  set) 
(change  from  hard  to  soft  contact  for  numerical  reasons) 
(contact  begins  when  gap  <  0.0625”;  contact 

pressure=5000psi  @  zero  gap) 

(specify  contact  surface  is  initially  bonded  and  may 
debond) 

(defines  rigid  surface  with  which 


elements  in  SYMMREST  are  to  interact;  this  surface  is 
coincident  with 

center  line,  extending  from  -1.0  to  13.0  inches) 

(names  node  1  a  set  of  nodes  called  CRACKREF) 

(a  crack  will  grow  ahead  of  the  plow  tip,  beginning  at 
node  1) 

(specifies  boundary  conditions  on  various  sets  of  node 
points) 

(restrains  y-displacement  at  nodes  along  stop  end) 
(restrains  x-displacement  at  nodes  along  long  edge  of  soil 
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PLOWNODE,  1 

PLOWNODE, 3, 6 
SYMMREST,1,6 

*RESTART,WRITE,FREQ=50 

*  STEP,NLGEOM,INC=  1000 

PLOW  ENTERING  SOIL  AT  5  FPS 

*  SOILS  ,CONSOLIDATION,UTOL=100 
0.0002, 0.1 6,  l.E-6, 0.01 


*BOUNDARY,TYPE=VELOCITY 

PLOWNODE,2,2,60. 

♦BOUNDARY 
EDGE, 8,8,0. 

♦CONTROLS, PARAMETERS=FIELD,FIEL 

D=DISPLACEMENT 

0.01,0.5 

*PRINT,CONTACT=YES 

♦DEBOND, ELSET=SYMMREST 
0.0,1.0,0.00067,0.0 

♦CRACK 

GROWTH, NSET=CRACKREF, TIME 
TOLERANCE=0.00067 
0.0,0.0,1.0,60.1 
♦NODE 

FILE,NSET=PLOWNODE,FREQ=l  0 
U,RF 

♦NODE  FILE,NSET=CLPOR,FREQ=l  0 

POR 

♦EL 

FILE,ELSET=CLSTR,POSITION=CENTRO 

IDAL,FREQ=10 

S,E,SINV,PE 

EVOL 

♦MONITOR, NODE= 1 000 1  ,DOF=2 


region) 

(restrains  rigid  body  displacement  of  plow  in 
the  x-direction,  z-direction,  and  about 
the  x-,  y-  and  z-axes) 

(completely  restrains  rigid  body  displacement  of 
centerline  plow  guide) 

(save  and  reuse  data  and  analysis  results  every  50 
increments) 

(begin  history  section;  large  displ.  analysis;  1000 
increments  allowed) 

(step  title) 

(effective  stress  analysis  for  fluid-filled  porous  media; 
transient  consolidation  analysis;  max  pore  pressure 
change  in  any  incr.  is  100  psi;  initial  time  incr.;  duration 
of  analysis;  min  and  max  incr.  sizes) 

(prescribe  constant  velocity  boundary  condition  for 
plownode;  in  y-direction;  60  in/s) 

(prescribe  zero  pore  pressure  boundary  condition  for 
edge  nodes) 

(relax  convergence  criteria  on 

displacement  field  equilibrium;  defaults  are: 
residual-0. 005;  solution  correction=0.01) 

(request  report  on  contact  states  in  the  message  output 
file) 

(specify  debond  and  crack  length  vs.  time  curves) 

(liner  decay  in  debond  force  with  time,  vanishing  after 
0,00067  s) 

(crack  length  increases 


linearly  to  60.1  inches  in  1  s,  beginning  at  node  crackref) 
(define  results  file  output  data  for  PLOWNODE) 

(displacement  and  reaction  forces  printed  every  10 
increments) 

(defines  results  file  output  data  for  node  set  CLP  OR) 

(pore  pressure  printed  every  10  increments) 

(define  results  file  output  data  for 


ele  set  CLSTR;  stress  and  strain  components;  stress  in 
variants;  plastic  strains; 
and  element  volumes) 

(monitor  plow  node  y-displacement  status  file) 
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*EL  PRINT, F=0 


(define  data  file  print  requests  for  element  variables; 
suppress  output) 

*NODE  PRINT, F=0  (define  data  file  print  requests  for  nodal  variables; 

suppress  output) 

*END  STEP  (end  of  step  and  analysis) 
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